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I. INTRODUCTION 



Numerical methods have been useful in the prediction of 
atmospheric flow patterns since the pioneering work of 
Charney et al. (1950) on the ENIAC computer. In that first 
experiment, and in operational models that soon followed, 
pressure height analyses were used as initial conditions for 
the prediction equations. At that time, the numerical 
models excluded gravity waves, and no special processing of 
the initial conditions were necessary. However, when the 
primitive equation (PE) models came into widespread use, the 
independent variables consisted of mass (pressure, tempera- 
ture and geopotential height) and motion (wind, vorticity 
and divergence) variables. Rather than analyze separately 
the wind field, it was computed from the geopotential height 
analyses with a diagnostic field relationship such as the 
balance equation (see Charney, 1955). Since PE models 
permit gravity waves as well as Rossby waves, unbalanced 
initial conditions may become badly distorted during the 
integration. Such imbalances occur when the motion and mass 
variables are not dynamically matched, as is typically the 
case due to the inaccuracies and the spatial distribution of 
the observations. Therefore, although winds are analyzed 
along with the pressure heights, some method is required to 
remove the gravity wave noise. 



With the advent of large quantities of asynoptic data 
(data observed at random times, such as from satellites and 
aircraft) came the concept of data assimilation, in which 
the numerical model plays a crucial role of carrying infor- 
mation between observation times (Charney, Halem and 
Jastrow, 1969; Hayden, 1973; Bengtsson and Gustavsson, 1971, 
1972; and Williamson and Kasahara, 1971). The primary 
motivation for data assimilation has been to update the 
numerical model frequently enough that it would represent 
the state of the atmosphere at all times. In this way, the 
asynoptic data could be used in the analysis more effec- 
tively, and thereby improve forecasts. However, several 
difficult problems became apparent from the initial, some- 
what naive, methods of updating a model. The most difficult 
problem has been the inability to assimilate all types of 
data. The mass observations are particularly difficult 
(Kistler and McPherson, 1975; Daley and Puri, 1980; Puri, 
1981) because the geostrophic adjustment mechanism of the 
model tends to disperse unbalanced mass information through 
the gravity waves. Another problem, particularly noticeable 
in the tropics, is that the model may become severely unbal- 
anced. This may set up spurious oscillations on a global 
scale that are very difficult to damp with time filters. 

For example, the pressure field in the tropics may tend to 
slosh with periods of 12 to 18 hours and amplitudes in 
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excess of 10 mb. It then becomes extremely difficult to get 
the model state to approach asymtotical ly that of the 
atmosphere . 

Balancing the analyzed data, which removes the gravity 
wave noise, avoids the problems mentioned above. Because 
the wind data may be located in regions without mass data 
and vice versa, the balancing method should be flexible. 
Although much of the gravity wave noise that interferes with 
data assimilation can be controlled with some sort of filter 
applied during the integration, most operational centers 
constrain the initial data so that winds and mass are 
bal anced . 

The two basic categories of the balancing procedure are 
the dynamic and the static. The dynamic method involves 
continuously inserting data until the model fields are 
adjusted to the new information. Integration may be 
performed in a forward-backward fashion as discussed by 
Nitta and Hovermale (1969), Temperton (1976), and Haltiner 
and McCollough (1975), or it may be performed in only one 
direction as discussed by Miyakoda et al. (1976) and Hoke 
and Anthes (1976). The main difficulty with the dynamic 
method is inefficiency; it requires the equivalent of a 24- 
hour forecast, and it does not seem to work well for mass 
data (Williamson and Temperton, 1981). Static methods, on 
the other hand, are widely used in operational centers 



(Daley, 1981). In the static method, a diagnostic relation- 
ship is imposed on the analyzed heights and winds. The 
acceptance of the mass data in the model depends on the 
constraints imposed during the balancing (Daley, 1978) or 
whether the analysis of mass variables also produces correc- 
tions to the motion variables (Philips, 1982b). 

Several static initialization methods are available. 

One method utilizes the multivariate optimum interpolation 
to make mass corrections consistent with motion corrections 
(Lorenc, 1981; Schlatter, 1975), and then the remaining 
gravity noise is removed with some sort of balancing. 
Unfortunately, this multivariate analysis method links the 
mass and motion through the geostrophic approximation, and 
therefore produces a bias around wel 1 -developed systems 
(Williamson, Daley and Schlatter, 1981). Another method 
requires that the mass and motion variables be analyzed 
independently, and then the calculus of variations initiali- 
zation of Sasaki (1958) either adjusts the mass variables to 
the motion variables, or vice versa, depending on the 
expected accuracy in the mass variables relative to the 
winds. The main difficulty with this method is that there 
is no convergence guarantee for the iterative methods 
required to solve these problems (Tribbia, 1981; 1982). 
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The constraints imposed on the initial data to remove 
gravity waves are most commonly the nonlinear normal mode 
methods of Machenhauer (1977) and Baer and Tribbia (1977). 

In these methods, the nonlinear component of the balance is 
computed assuming little or no tendency of the gravity mode 
coefficients. No other initialization method is capable of 
suppressing the initial imbalance in a forecast so effec- 
tively. Additionally, the nonlinear normal mode methods 
have the advantages that they provide conditions that are 
compatible with the numerical scheme of the model, generate 
realistic vertical motion in the extratropics, and produce 
balanced flow in the regions with terrain (Daley, 1981). 

In a theoretical study, Leith (1980) used a 
quasigeostrophic model to show that the nonlinear normal 
mode methods are nearly the same as constraining the initial 
conditions to the balance and omega equations. Therefore, 
it seems reasonable to expect that the balance equation 
might still be competitive with the normal mode methods, 
particularly since the balance equation constraint is 
relatively easy to impose. 

There remains the problem of generating an appropriate 
divergence to go with the nondivergent winds produced by the 
balance equation. Tarbell et al. (1981) used a modified 
omega equation which improved the precipitation forecasts of 
a mesoscale model during the initial hours. Considering 
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that the omega equation, and even the nonlinear normal mode 
method do not generate divergence patterns in the tropics 
that are compatible with the latent heat release (Bengtsson, 
1981), the divergence from the forecast first-guess may be 
the best estimate for divergence. 

In this study, we examine the adequacy of the forecast 
first-guess divergence. Because this divergence is gene- 
rated by the model, it is compatible. Thus, the interruption 
to the cumulus convection in the tropics is minimized during 
assimilation. Another goal of this study is to determine 
whether the classical balance equation method of Charney 
(1955) used in the variational method of Sasaki (1958) is 
competitive with the more elegant, yet cumbersome, nonlinear 
normal mode method. Since one of the main benefits of the 
nonlinear normal mode method is the divergence it generates, 
comparison of the balance equation and normal mode method 
is, in many respects, a test of the accuracy of the forecast 
first-guess divergence. 

In the following, a brief description of the data 
assimilation system used in this study is given. Besides 
the initialization methods, this system includes an 
objective analysis method for wind and pressure height 
observations and the global finite difference model 
described by Arakawa and Lamb (1977). The initialization 
methods, which are the primary topics of this report, are 
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presented in detail. Both methods are presented in a 
calculus of variations framework. The normal mode method 
uses Machenhauer 's (1977) initialization as a constraint, 
whereas the other method uses the balance equation as a 
constraint. Finally, several data assimilation experiments 
are described that illustrate some of the characteristics of 
the different initialization methods, particularly with 
regard to precipitation rates and divergence during the 
early forecast hours. 
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II. THE DATA ASSIMILATION SYSTEM 



In the data assimilation method used in this study, the 
prediction model is periodically updated at 1 2 - h intervals. 
Each update requires several steps. First, the 12-h fore- 
cast is interpolated to the analysis coordinates, which are 
pressure surfaces on a 2.5° by 2.5° grid. This forecast 
becomes the first-guess field. The objective analyses of 
wind and pressure height are done with a three-dimensional 
successive corrections method (See Appendix A) on the 
standard pressure surfaces (50, 100, 150, 200, 250, 300, 

400, 500, 700, 850 and 925 mb). However, the surface 
pressure analysis is produced from the calculus of varia- 
tions method of Holl and Mendenhall (1972). At this time, 

% 

the initialization is conducted. The balance equation 
initialization is done on pressure surfaces prior to the 
interpolation to model coordinates, whereas the normal mode 
initialization is done in model coordinates. These initial- 
ization methods are described in the next two chapters. 
Finally, a 12-h forecast is made in preparation for the next 
update . 

The other variables in the model, such as boundary 
layer depth and moisture, are not updated with new data. 
Rather, they are carried forward from the forecast first 
guess. 
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The numerical model used to assimilate the data was 
developed at the University of California at Los Angeles 
(UCLA) to study the general circulation of the atmosphere 
(Arakawa and Mintz, 1975; Arakawa and Lamb, 1977). In this 
model, the resolution of the mass variables (surface pressure 
and temperature) is 4° latitude by 5° longitude and six 
levels from 50 mb to the surface. The zonal wind is defined 
one-half grid interval to the east and the meridional wind 
is defined one-half grid interval to the south (Scheme C in 
Arakawa and Mintz, 1975). The heating parameterization 
includes the Arakawa and Schubert (1974) parameterization 
interacting with a bulk parameter boundary layer (Randall, 
1976; Lord, 1978). 

The time differencing is a combination of five leap- 
frog steps for each Matsuno backward step, while the heating 
is computed during a single forward step preceding each 
Matsuno step. 

In the next section, the model's normal modes are 
computed, and more detail is given on the dynamical forecast 
scheme . 
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III. NONLINEAR NORMAL MODE INITIALIZATION 



Before nonlinear normal mode initialization can be 
applied, the normal modes of the linearized model equations 
are required. These are obtained below by separation of 
variables of the linearized equations. Unfortunately, the 
normal mode methods match the motion variables more closely 
than the mass variables (Daley, 1981). To establish control 
over this mass rejection mechanism, the normal mode method is 
converted to a variational one similar to that of Daley 
(1978) and Tribbia (1982) . 



A. VERTICAL MODES 

The linearized governing equations may be written with 
the vertical component in vector form as 



\v + f ikx \v + 


7(RT In p s + $) = Q w . 


(3.1 ) 


T + t (v* V) 


= Q t’ 


(3.2) 




- 1 




In p $ + n T ( 7 . 


V) = Q p 5 and 


(3.3) 



<|) - $ + G T 

Here W is the vector form of wind,T is temperature and T is 
the rest state temperature. P s is surface pressure, <f, is 
geopotential, v • IV i s divergence and <j) S is terrain 
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geopotential. t, H and G are linearized operators defined 
in Appendix 3. Q , Qj and Qp are the nonlinear components 
of their respective equations. 

Following Temperton and Williamson (1981), the 
temperature and surface pressure may be described by a 
single variable by using 

gh = <p + RT In p (3.5) 

Operating on (3.2) with G, multiplying (3.3) by RI, and then 
adding the resulting two equations gives a single equation 
for mass, 

g ^ h + G x(v : W) + RT T ] ( 7 • W) = G Q T + Q p RT ( 3 . 6 ) 
which may be rewritten as 

g TY h + C (7- \V) = g Q h ( 3 - 7 ) 



where 

C = G - + R T n T (3.8) 

and 

g h ■ G Q t + R Q p T (3.9) 

The equation set (3.7) is vertically coupled, but by 
separation of variables, it can be transformed into a set 
that is not coupled. This is done through the identity 

E' 1 C E = g D (3.10) 
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where matrix i and diagonal matrix 0 contain the eigen- 
vectors and eigenvalues of C, respectively. Transforming h 



a n d \V , or 

W = E' 1 \V and (3.11) 

h = E _1 h , (3.12) 

produces the following uncoupled equation set: 

fr V + flk • 7X \V + g (h) = Q w and (3.13) 

ft h + 0 (7- \V) = Q h (3.14) 



where and Q h are transforms of and Q h , respectively. 

The independent variables in (3.13) and (3.14) are the 
coefficients of the vertical modes., i.e., the eigenvectors 
£. Naturally, there are as many modes as there are levels 
in the model. 

The eigenvectors, E, are shown in Figures 1 and 2 for T 
equal to a constant (300°K) and for T equal to (209, 214, 
233, 254, 270, 2 3 1 ) 0 K . The corresponding eigenvalues 
(equivalent depths) are given in Table 1. Notice that while 
the profiles do not significantly change shape for the 
different temperature profiles, their eigenvalues are 
considerably different. Also, the vertical modes of the 
Arakawa and Lamb (1977) model shown here are noticeably 
different from those given by Temperton and Williamson 
(1981) for the same values of T. Specifically, the peaks in 
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Fig. 2. Same as Fig. 1 except for a rest-state temperature of (209, 214, 
233, 254, 270, 281)°K. 



Table 1. Equivalent depths (eigenvalues) of the vertical 
modes of the Arakawa and Lamb (1977) model. 



Level 


T Constant 
300°K 


T Average of 
1 March 1965 


1 


9,291 m 


7,874 m 


2 


1,689 


784 


3 


387 


171 


4 


126 


55 


5 


45 


18 


6 


10 


3 



the profiles near the model top are not present in the 
gravest modes given in Figures 1 and 2. Since these modes 
are insensitive to the values of rest state temperature, 
they need not be updated for each new data set. 

The equivalent depths are sensitive to the location of 
the model top. For example, changing the model top from 50 
mb to 0 mb increases the equivalent depths of the external 
mode from 7874 m to 9660 m. The consequence of keeping the 
model top at 50 mb is that all the vertical eigenvalues or 
equivalent depths are smaller than if the top was placed at, 
say, 10 or 20 mb. This becomes a factor when the nonlinear 
balance is performed, as will be discussed later. 

After h is balanced, the inverse of (3.5) is required 
to solve for P $ and T. This is done following the approach 
of Temperton and Williamson (1981) and Andersen (1977), 
where v*\V may be eliminated between (3.3) and (3.7) to give 
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The nonlinear components are grouped in the term NLT. From 
linear theory, it is possible to relate h adjustments to 
p s adjustments, or 

A In p s = g Ah), (3.16) 

where the adjustment balances the analyzed fields. 

In a similar manner, elimination of v-\V between (3.2) 
and (3.7) produces 

AT = t [C" 1 ] g Ah (3.17) 

Using definition (3.5), it is easily shown that (3.16) and 
(3.17) are consistent methods for determining Alnp s and aT 
from Ah. 

A two-grid interval wave was found to exist in the 
solution of (3.17). To eliminate this problem, a varia- 
tional method was used wTiich minimizes aT, thereby removing 
the two-grid length waves. In this method, a<j> is computed 
from (3.5) and then aT is determined from the minimization 
of 



N 2 2 N 

I (AT L ) 6 ( A<f> L ) + 2Xy(A4)|^- E ^l^aT^Ac^. 

— 1 K — 1 



(3.18) 



The details of this solution are found in Barker (1981c). 
This method reduced the size of the root mean square (RMS) 
corrections to about one half of those obtained using 



In this section, the vertical coupling of the 
linearized equations (3.1)— (3.4) was eliminated through 
separation of variables. Following a similar procedure, the 
horizontal coupling can also be removed, as shown in the 
next section. 



8. HORIZONTAL MODES 

The equations (3.13) and (3.14) in linearized finite 
difference form for each vertical mode are 

, - , (3.19) 

~~Z ^ » 6 g(<5.h). + i /2 ■; 

• - f . (v) /0 . + * — i ZlLLil = Q 



6^U 



t i + 1 / 2 , j " 'j vv; i+l/2,j 
— 0 



a a - AX 

J 



u i + 1 / 2 , j , 



6 t v i,j-T/2 



+ (fg u) i,j-1/2 + g(lS 9 h) i , j + 1 / 2 = 

a j-l/2 aA0 



(3.20) 



i » j - 1 / 2 



and 



s t h i,j 



+ 0 



aA0 






a j a A X 



= Q 



h i , j 



(3.21) 



The mode index, m, is assumed for D and each variable. 
The finite difference operators are 



( 6 x T )k = T k + l/2 + T k -1/2 » 



and 




T k+1 /2 + T k-1/2 
2 




(T) k • 
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The other variables are defined as follows: a is the earth's 
radius, AA is the longitudinal grid interval, A0 is latitu- 
dinal grid interval, u and v are the east and north compo- 
nents of wind, respectively, i and j are the longitudinal 
and latitudinal indexes, respectively, and cr is cose. 

The linearized equations ( 3.1 9) -( 3.21) are derived from 
the model equations in the flux form given by Arakawa and 
Lamb (1977). However, linearization of these equations 
removes their flux form, and the only terms absolutely 
unique to the model are the Coriolis terms, which differ 
from other models in the way that f is averaged. As it 
turns out, special definitions for the Coriolis terms are 
desirable in order to keep the matrix operator of (3.19)- 
(3.21) symmetric, which simplifies the corresponding 
eigenvector matrix so that it can be inverted with a 
transpose operation. This decreases the computer storage 
and time requirements needed to transpose the variables into 
normal mode coefficients to one half that required for non- 
symmetric operators. Fortuitously, the small errors 
introduced by these modifications do not affect the normal 
mode balancing (Temperton and Williamson, 1979). To achieve 
symmetry, the Coriolis term in (3.19) is replaced by 

— — A + — t“A 

f i a ,i-1/2 (v) i+l/2,j-1/2 + f ,j J j+l/2 (v) i+l/2,j+l/2 

2 a ■ 
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and the Coriolis term in (_3.2Q) is replaced by 





2 



where 



f . 
J 




and 




These definitions correspond to a potential enstrophy 
conserving finite difference scheme as derived by Temperton 
and Williamson (1979). 

From this point, the procedure closely follows that of 
Dickenson and Williamson (1972) except that the finite 
differences are written for scheme C (see Arakawa and Lamb, 
1977; Temperton and Williamson, 1981). The dynamical state 
vector is defined as 




(3.22) 



By assuming a wave solution in the form 



1-1 . 

y(x. , e • ,m ) = z y ( k , 9 • ,m)e 

- i J k=0 ~ J 



i kx • 



(3.23) 
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and consequently, 



yC^, 0 . ,m) = -=- z y ,e. ,m)e 
J 1 i=l - 1 J 

(3. 1 9) -(3. 21 ) become 



-ikx 






(3.24) 



A u.m'(k) - ^ [ f j <J j-i/2 v j-l/2 + f j a j + 1 / 2 v i + l/2 ] 



i k 1 ^ 

g— h. s, = Q , 
3 aa. j f x u * 

J 



(3.25) 



IT ^ j -l / 2 + ^ [ f j : ij + fj-i ij.ia-’tk) 



(3.26) 



aA9 ^ h j“ h j-l^ = Q v and 
D - 1 



3t 



j + aa. tlk ' U j S f + A0 (v j+l/2 a j-l/2 

V j-l/2 a j-l/2 )] = V 



(3.27) 



where 

m'(k) = cos(^-)-i sin(^-), 

k ' = sin(^)/(^-) , r(k) = cos(^-), and S f is the 

filtering applied near the poles to keep the gravity-wave 
terms computationally stable (Arakawa and Lamb, 1977). 

The symmetric operator 

^m 1 ( k ) 0 0 



m 



0 

0 



-l 



0 



0 ( 9 /D ra )1/2 / 
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is used to redefine 

Y » 1 .e. , Y 1 = S m Y 



( 3 . 28 ) 



so that 



- io j JT'V ' ^ tf ?j-l/2 } j-t/2 * f jVl/2Vl/2 ) + 



f k h J S f 1 % 



I ~ I - I 



( 3 . 29 ) 



1 0 j -1 / 2 i7 (v j-l/2 ) " 0 0 - 1 / 2 



£ f j * f j + -i a i-i ] • 



°j -i /2 riV [h j ' h j-i ] = 



( 3 . 30 ) 



• ic j sT (i V + T 1 k “j s f + 546 [v j+l /2 Vl/2 



V • 



j -1/2 j - 1 / 2 



] = Q, 



( 3 . 31 ) 



or in matrix form, 



1 9 sT i’ + M 



- iQ H ' 



( 3 . 32 ) 
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0 



Q = 



where 



•Q, 



o 

0 



L = 



a j -1 / 2 






0 



2 

0 A 



and 



N B N 



For the general case when j is not near the north pole or 
the equator 



a j ■ 



r ( k ) 

2" f j ° j - 1 / 2 



_ a r lil ) f " 

a j -1 / 2 2 r j 



L 



J!L k S 
a K b f 



m 



aA9 j -1/2 



+ J» u c 
a K b f 



m 



j - 1 / 2 aA9 
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B 

J 



and 



0 

0 



. rlLL f + 

2 J 
0 



Vl/2 



m 

aA9 °j+l/2 



0 

0 

0 



0 



_a j - 1 / 2 



rm. 

2 




0 0 ‘ 

C 

n „ _m_ 

0 a j-l/2 aA9 



0 



0 



0 



To produce continuous Coriolis forcing near the poles, 
a computational u is defined. At the North Pole, 

[u] N = 0, where 



the [ ] represents a zonal average. 
Continuity requires that 



(u i+l/2,N 



, Sf a N - 1 / 2 

U i-1/2,N ) aX " V i,N-l/2 (A9/2) = 



rvl °N - 1 / 2 

L J N-l/2 ( a 9 / 2 ) 



(3.33) 



It is important to notice that except for wave number zero 
(MO), 

Mn-1/2 = 0- 



Consequently, L is one row and column larger when K is zero 
than it is otherwise. Transforming the variables using 
(3.23) gives 



U N “ av N-l/2 



(3.34) 
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where ot — 2 j ^ 1/2 ^ k a 0 ) • 

The v - equation becomes 

^ * r(k) _ - i 

1 CT N - 1 / 2 Jt v N - 1 / 2 2 a n -1 / 2 ^ f N aV N-l/2 + 

(3.35) 

f + '■ n g N- 1 / 2 r r r' r ' _ z ’ 

f N- 1 U N-1 ] ' aAG C m*- h N " h N-l ^ " Q v * 

The continuity equation at the pole is 

^ h M z' 

W D n/% ds/ flrea ’ 
where Area = I(aA0 a i^_i /2^ 4 

and 

r „ i 

? v N ds = “N aAe c N-l/2 • 



h^ is zero except for wave number zero, so that 

I 

i 3 - 1 ^N-l/2 ~ 1 

" 4 CT N - 1 / 2 9 t h N " C m “aAG a N-l/2 " Q h * 



Therefore (3.35) and (3.36) give 



A, 



C m a N-l/2 



0 aAG 

C m CT N-l/2 0 

aAG 

r(k) f + 

2 CT N - 1 / 2 r N-l 

0 



N-l /2 c 
aAG m 

0 



and 



B 



N-l 




(3.36) 
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A grid with a four degree meridional interval puts the 
equator at a v-point. Using v-j^ as th e equator point, the 
equations can be written for the symmetric and antisymmetric 
components. For the symmetric component, u-j = u Q , h-| = h Q 
and v-j/2 = The h- and u-equations are 



i _L r( k) r -+ 

" 1 a l at U 1 2 *- f i a 



1+1/2 V l+1/2 ] + 



-f k ‘*l S f = ^u and 

' 1 a l il V k u . S f + 5 a? a l +1 / 2^1 +1 / 2 



( 3 . 37 ) 



= Q 



h ’ 

( 3 . 38 ) 



so that 



’0 

Ea k ' s 

_ a k b f 
0 

0 



-J5- k S 






T~ f l a l+l/2 



a l +1 / 2 



m 

aAe 



-i o-i 

o o,J 



and C2 = B-j"*" . 

For the antisymmetric component, u-| = -Ug, h-| = -hg and v-^2 
is not identically zero. The equations are 
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1 a l 3 1 U 1 ' {f l a 1 / 2 ^1/2 + f l a l+l/2^1+l/2 } + 



-f k fi l S f - ’ 



(3.39) 



1 a l/2 aT ^1/2 ' a l / 2 {f l U 1 " f l ( " u l )} 

a 1 / 2 iZe ( *1 + M = Si and 



(3.40) 



1 a lTt ^1 + a k U 1 S f + 



(3.41) 



lA? { Vl/2 a l+l/2 " V l/2 a l/2 } = Q h ’ 



o r 



_ r(k)_ f - 

a l / 2 2 f l 



-r- k S 



r ( k ) 

“T” f l a l / 2 



- CT 



m 



1/2 aA0 



3 k s 
a K b f 



a l / 2 
' aA0 



0 

0 



r( k ) f + - 

r l a l+l/2 



m 

aA0 a l +1 / 2 



0 

0 



and 
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0 



0 



Qi = 



°l/2 /2 



0 



Changing the model resolution will generally require that 
these equations also be changed, depending on how the 
variables are staggered relative to the equator. 

To complete the computation of the horizontal modes, 
the matrix equation (3.32) is rescaled using 
y = Q^ 2 y', H = Q^/ 2 H ' 

and 



L - Q- 1 / 2 L Q- 1 / 2 . 

This allows (3.32) to be written 

- i gi/2 £ i ♦ qi/2 L j < - 




or 



+ = (3,42) 
If Y contains the eigenvectors or normal modes of L, then 

this equation becomes 

- i J_ [ Y " 1 y ] + Y' 1 L Y ( Y -1 y ) = - i Y _1 H. (3.43) 

3 1 * ~ ~ ~ — 

The identity 

A = Y’ 1 L Y , (3.44) 

where A is a diagonal matrix containing the eigenvalues of 
L, makes it possible to rewrite (3.43) as 



J-C = -\liC + r 

n "■ - 



(3.45) 
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This is the wave equation of the expansion coefficients, C, 
that we have been seeking. The elements of C are functions 
of the vertical, zonal and meridional mode numbers m, k and 
1, respectively. These coefficients are the amplitudes of 
the various modes required to represent a particular atmos- 
pheric state, i .e., 

C = Y' 1 y . (3.46) 

The nonlinear term is now r, and the mode frequencies are A. 
For each m and k, there are 3N equations for C; 2N are 
gravity waves and N are rotational waves. N, in this case, 
is the number of degrees of freedom in the meridional 
direction . 

The structure of the modes from (3.43) is given in 
Figures .3 and 4. They are nearly identical to the modes 
published by Dickenson and Williamson (1972). The scaling 
is different, but this is of no consequence, as the modes 
are put in orthonormal form before they are used. The 
frequencies of the various modes corresponding to those 
given by Dickenson and Williamson (1972) and Temperton and 
Williamson (1981) are given in Tables 2 and 3. Resolution 
differences account for most of the variability of the 
frequencies, which can be seen from the computations of 
Dickenson and Williamson (1972) for two different resolu- 
tions. Their 5° unstaggered grid results are similar to the 



42 



0.5 



0.5 




k=A, 




0=10 km 




.0.3 H II I I M I IN I M 1 1 1 

0 20 40 60 80 

Latitude 



0.5 



k=l , 1 R =3 , 0=1 0 km 




-0.5 



JJJJ-LLLLLi Ml l 1 ll i 



0 20 40 sa so 

Latitude 




Latitude 



Fig. 3. Structure of selected rotational modes for the 
model used in this study, which may be compared 
directly with those of Dickenson and Williamson 
(1972). _ 
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Fig. 4. Similar to Fig. 3 except for selected 
gravitational modes. 
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2. Frequencies of Rossby modes for D = 10 km, k=l 
from computations of Temperton and Williamson 
(1981) (T&W), Dickenson and Williamson (1972) 
(D&W), and for the model used in this study, B. 
The horizontal grid intervals are specified in 
degrees . 



T&W , 


10° 


D&W, 


5° 


D&W, 


2.5° 


B 4°x5° 


6.11 


E-05 


6.12 


E-05 


6.14 


E-05 


6.14 


E-05 


1.44 


E-05 


1.43 


E-05 


1.45 


E-05 


1.45 


E-05 


8.64 


E-06 


8.60 


E-06 


8.73 


E-06 


8.74 


E-06 


5.72 


E-06 


5.73 


E-06 


5.87 


E-06 


5.88 


E-06 


3.98 


E-06 


4.02 


E-06 


4.17 


E-06 


4.17 


E-06 


2.87 


E-06 


2.93 


E-06 


3.08 


E-06 


3.09 


E-06 


2.14 


E-06 


2.20 


E-06 


2.36 


E-06 


2.36 


E-06 


1.63 


E-06 


1.69 


E-06 


1.86 


E-06 


1.86 


E-06 


1.27 


E-06 


1.32 


E-06 


1.49 


E-06 


1.49 


E-06 


1.01 


E-06 


1.04 


E-06 


1.22 


E-06 


1.22 


E-06 


8.10 


E-07 


3.18 


E-06 


1.02 


E-06 


1.02 


E-06 


6.62 


E — 0 7 


6.43 


E-06 


8.58 


E-07 


8.57 


E-07 


5.52 


E-07 


4.99 


E-07 


7.30 


E-07 


7.30 


E-07 


4.70 


E-07 


3.77 


E-07 


6.27 


E-07 


6.28 


E-07 


4.11 


E-07 


2.71 


E-07 


5.43 


E-07 


5.45 


E-07 


3.75 


E-07 


1.75 


E-07 


4.73 


E-07 


4.76 


E-07 


3.13 


E-07 


8.60 


E -08 


4.14 


E-07 


4.18 


E-07 
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Table 3. Similar to Table 2 except for frequencies of 







eastward 


gravity modes 


for D 


= 10 


km, k = 1 . 




1 


T&W , 


10° 


D&W , 


5° 


D&W, 


2.5° 


B 4°x5° 


0 


-5.44 


E-05 


-5.38 


E-05 


-5.38 


E-05 


-5.38 


E-05 


1 


-1.31 


E-04 


-1.30 


E-04 


-1.30 


E-04 


-1.30 


E-04 


2 


-1.87 


E-04 


-1.85 


E-04 


-1.86 


E-04 


-1.87 


E-04 


3 


-2.35 


E-04 


-2.33 


E-04 


-2.36 


E-04 


-2.36 


E-04 


4 


-2.79 


E-04 


-2.78 


E-04 


-2.83 


E-04 


-2.84 


E-04 


5 


-3.22 


E-04 


-3.20 


E-04 


-3.29 


E-04 


-3.31 


E-04 


6 


-3.63 


E-04 


-3.61 


E-04 


-3.75 


E-04 


-3.78 


E-04 


7 


-4.01 


E-04 


-4.00 


E-04 


-4.21 


E-04 


-4.25 


E-04 


8 


-4.36 


E-04 


-4.34 


E-04 


-4.66 


E-04 


-4.71 


E-04 


9 


-4.69 


E-04 


-4.67 


E-04 


-5.10 


E-04 


-5.18 


E-04 


10 


-4.98 


E-04 


-4.96 


E-04 


-5.54 


E-04 


-5.64 


E-04 


11 


-5.23 


E-04 


-5.21 


E-04 


-5.97 


E-04 


-6.09 


E-04 


12 


-5.44 


E-04 


-5.42 


E-04 


-6.38 


E-04 


-6.54 


E-04 


13 


-5.61 


E-04 


-5.59 


E-04 


-6.79 


E-04 


-6.99 


E-04 


14 


-5.72 


E-04 


-5.68 


E-04 


-7.19 


E-04 


-7.43 


E-04 


15 


-5.94 


E-04 


-5.75 


E-04 


-7.57 


E-04 


-7.86 


E-04 


16 


-5.94 


E-04 


-5.91 


E-04 


-7.94 


E-04 


-8.29 


E-04 
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Temperton and Williamson (1981) 10° staggered grid values, 
whereas their 2.5° nonstaggered grid results are similar to 
those computed from the above equations for a 4°X5° 
staggered grid. As Temperton and Williamson (1981) point 
out, the staggered grid produces modes comparable to those 
of a nonstaggered grid with half the resolution. 



C. NONLINEAR BALANCING 

Machenhauer (1977) discovered that the nonlinear terms 
have a much slower time variation than their respective 
gravity modes. Under this first order approximation, the 
equation for the gravity modes from (3.45), 

y|-C(k,A,m) = -ivC(k,A,m) + r(k,A,m) (3.47) 

has the solution /, * 0 ^ 



C( k , a , m, t ) 



r( k , A ,m) + 

iv(k ,A ,m) 



[C( k , A ,m , <j> ) 



r(k,A,m) -i-i vt 
i v J 



Therefore, removal of the fast modes requires that the 
second term on the right-hand side of (3.48) be zero, i.e.. 



C B (k,A ,m,i> ) = r ( k , A ,m) /( 1v ) . 



(3.49) 



The subscript 3 signifies the balance condition. 
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Rather than define the nonlinear terms, it is easier to 
use the model to make this computation. This is done by 
running the model for one time step and determining the 
tendency of the gravity mode coefficients. The nonlinear 
term is the total minus the linear tendency, or 



( k , £,m) = , m >°) + iv„C(k,a,m,t) . 



(3.50) 



v c is the computational frequency for the forward timestep, 
which may be determined from the analytic frequency by 



v = arctan 



(3.51 ) 



using standard methods. Using v c for v in (3.49) and then 
substituting (3.50) gives the corrections to the fast modes 
required to balance, or 



ACg = -i ( C ( k , i ,m , At ) -C ( k , i ,m ,0 ) ) / ( Atv c ( k , i ,m) ) (3.52) 

Leith (1981) used the quas i -geostroph i c relations to show 
that this balance is equivalent to using a slightly modified 
version of the balance equation along with the omega equa- 
tion. Phillips (1981) points out, however, that more than 
one iteration with (3.52) introduces new terms not 
consistent with quasi-geostrophic balancing, and therefore. 
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should be avoided. The Baer-Tribbia (1977) method, however, 
does allow more iterations, and will be tested in future 
versions of the normal mode balance. 

Additionally, Machenhauer's method (3.52) does not work 
except for those gravity modes with periods equal to or 
less than any of the rotational modes (Ballish, 1981; 
Phillips, 1981). For the version of model used here, this 
restriction limited the modes that could be balanced to all 
of the external and first internal, and a few of the second 
internal, modes . 

D. VARIATIONAL BALANCING 

One of the difficulties with the nonlinear normal mode 
method is that the rotational modes are primarily determined 
from the vorticity of the analyzed wind fields (Daley, 1981; 
Daley and Puri, 1981). This becomes a problem over vast 
regions of the globe in which the principal observations are 
remotely sensed temperature profiles from satellites. 
Therefore, to describe these regions adequately, the 
initialization must be able to assimilate mass observations. 
A possible method for incorporating mass observations is to 
convert the normal mode balance into a variational framework 
(Daley, 1978). Tribbia's (1982) spectral shallow water 
approach was adapted for this purpose, except here the 
method is developed for a multilevel finite difference 
model . 
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The normal modes introduced by (3.43) may be written 






U(i ,0j. ) 

Ul.Qj) 



exp(-ikx i )/I 1//2 



(3.53) 



but instead of (3.46), the coefficients of the rotational 
modes are determined from a specially designed inner product 

C(k ’ £) =<rr (k ’ £) > (3.54) 

where 

(3.55) 



<y-Y(k,i)> = 



I 

z 

i-1 



^ A i + 1 / 2 ^ u ^ 9 j ’^i + l/ 2^^ £,9 j^ a ^ 9 j ^ 

J ^ 



+ W u^ 9 j-l/2 ,X i^ V ^ e j-l/2’ A i^ V(A.0j)3a(0j. 1/2 ) 

- *1/2 
+ W^( 9 r X i )[^(9 j ,X i ) $ (£,9 j )]a(9 j )}exp(-ikX i )/I 

W,. and W, are the weights for winds and mass, respectively. 

9 

a is the same as in (3.25) — (3.27). Because no attempt is 
made to put a vertical variation in the weights, the 
vertical mode index has been dropped from these equations. 

The inverse of (3.46) is 

1-1 N „ 

Y = Z Z C £ Y £ (3.56) 

k=0 a»l K ~ K 
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If yj is the balanced initial condition and yo the 
analysis, then an optimal initial condition is that 
obtained from a minimization of 



1 = ' Io ),( Il - y 0 )> 



(3.57) 



with respect to the modal coefficients. 

The balanced data have components on both the rota- 
tional and gravitational manifolds, i.e., 

h = IR + iG’ (3-58) 

where 



IR = 



1-1 R 

Z Z 

k =0 4=1 




and 



(3.59) 



y G 



1-1 G 

z z 

k =0 4 = 1 



G k yf 
l 



(3.60) 



X k and G^ represent the rotational and gravitational mode 
coefficients, respectively. Notice that G plus R equals N, 
which is the total number of modes. The nonlinear balance 
relationship requires that the gravitational component be a 
function of the rotational component only (Phillips, 1981; 
3aer and Tribbia, 1977; Daley, 1981), i.e., G^ = G k (X), 
therefore minimizing (3.57) requires that 



il 

9 X 




a 



(3.61 ) 
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for the rotational mode coefficients X, or 



2 <( 





a 



(3.62) 



As Tribbia (1982) notes, this equation is not easily solved. 
This is particularly true for non-zero values of G. To make 
this equation more tractable, an iterative solution is 
possible where the first pass is solved using 

y[ ]) = (3.63) 

The superscript is iteration cycle. This simplifies (3.62) 
where 




(3.64) 



reduces to a linear equation set 



A X< ] > = Z 



(3.65) 




and 




for all Y s on the rotational manifold. 
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If the weights are unity everywhere, then the orthonormal i ty 
of the modes makes A an identity matrix, and yj is simply 
equal to the rotational component of yq. 

Subsequent values for X may be computed using (3.62). 

G is evaluated from Machenhauer' s equation (3.52) using the 
value of X from the previous cycle. The variation of G with 
respect to X s is computed numerically. The benefits of 

a 

following this complex procedure are almost certainly not 
worth the effort. 

To avoid the computational workload of (3.62), 
the variational balance is performed on the analyzed correc- 
tions. Assuming the corrections are much smaller than the 
total analysis, the approximate equation (3.64) is more 
accurate than it would be if the balance were performed on 
the total analysis. 

The solution is still not very easy. For total 
variability of the weights, the equation set (3.65) 
represents ( 1/2 • R ) or 1656 equations. The correspond i ng 
size of A is over 2.5 million elements. Therefore, trhe 
computation of the inverse of A is extremely cumbersome, 
especially when the computer memory available is only about 
100,000. Restricting the variations of the weight to 
latitude only decreases the size of A to 23X23. But now we 
must solve 216 different arrays, one for each vertical and 
horizontal mode configuration. These are solved only once 
and stored. 
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Some attempts were made to add longitudinal variation 
to the weights by iterating, but this did not work. 

Firstly, it was not possible to keep -yj on the rotational 
manifold. Secondly, the variations in the magnitude of the 
weight required to influence significantly the result were 
too large for convergence. 

Ultimately, it is hoped that a subset of the rotational 
manifold may be found that still defines the important 
meteorological information that needs to be retained by the 
initialization and yet has reasonable array sizes. 

A summary of how the final method works is given by a 
Leith (1981) manifold schematic in Fig. 5. Prior to the 
update, the model is assumed to be on the slow manifold (M), 
say at point 0. Adding the rotational component puts the 
model on the data manifold. This step is represented by the 
line between 0 and 1, and is equivalent to linearly 
balancing the corrections and inserting them into the model. 
Note that only the rotational modes are affected by this 
step, where (3.65) is solved for X and then the rotational 
component of the corrections is determined from (3.59). The 
imbalances introduced by updating are removed by solving 
(3.52) and then (3.60). This is the nonlinear step and is 
represented by the line between 1 and 2. 
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Fig. 5. Manifold schematic of the update and balance 

procedure used in the data assimilation. Point 0 
represents the forecast before the update. The 
path between 0 and 1 represents the linear balance 
step. The path between 1 and 2 represents the 
nonlinear balance step. 
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E. TESTS OF THE NORMAL MODE METHOO 

The magnitude of the rotational and gravitational com- 
ponents of the analyzed corrections were computed from real 
data taken from the 00 GMT 16 Nov 1979 analysis (see Figs. 6 
and 7). The surprising thing about these curves is that the 
gravitational component is as large or larger than the 
rotational component. This is especially true of the 
surface pressure. In fact, the surface pressure gravita- 
tional component is so large that the two components must be 
of opposite signs over much of the globe. Because the 
nonlinear component of the initialized fields is determined 
from the rotational component only, the gravitational 
component of the analysis is not used and therefore repre- 
sents lost information. It is also a measure of how well 
the analysis dynamically matches the wind to the mass. 

Completeness theorems for the normal modes require that 
the sum of the rotational and gravitational components of 
the data be equal to the original data. As a test of the 
approximations and final code, the degree to which the 
completeness theorem applied was determined. This was done 
by computing the rotational component of the corrections, 
and then computing the gravitational component from what 
remained of the corrections after the rotational component 
was removed. The corrections were then reconstructed by 
adding the rotational and gravitational components. These 
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Fig. 6. Update analysis (solid line) and its rotational component (dashed 
line). Subscript 3 designates level 3. 
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reconstructed data were compared to the original data by 
computing RMS differences (see Fig. 7 ) . The differences 
between the two wind variables were essentially zero, but 
nonlinear relationships between surface pressure and temper- 
ature caused significant residuals to remain in the mass 
variables. The long dark dashed lines near the abscissa in 
Fig. 7 (c) and (d) show the residual for level three of the 
model and the numbers above the curves give the global RMS 
averages. It can be seen that the temperature residual is 
largest, and represents about 14% of the total correction. 
The pressure residual is only about 2%, however. 

As a test of the variational technique, the linear 
balance was rerun with different weights on the mass varia- 
ble each time. Figure 8 shows the global RMS residuals 
between the balanced and analyzed variables for the differ- 
ent weight values. Clearly, the most sensitive region is 
for weight values between 0 and 1. For a loss in the wind 
residual of about 0.4 m sec"-*-, the improvement in the mass 
variables is 0.3°C and 1 mb. Increasing the weight from 1 
to 10 does not dramatically improve the fit of the balanced 
mass variables to the analysis. Eventually, for very large 
weight on geopotential, the wind residual becomes as poor as 
the forecast first guess. It is interesting to note the 
limit in the residual of the wind variables. Even when the 
weight on the mass variables is zero, the wind residual is 
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Fig. 8. The RMS differences between the analyzed and 

initialized variables for different weighting on 
the mass analysis, W<$. The averages are for the 
entire globe and all levels. The dashed straight 
line shows the RMS differences between the 
analysis and the forecast first-guess. 
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2.8 m sec' 1 . The reason for this is that the analyzed wind 
contains unrealistic divergence that is orthogonal to the 
rotational manifold (Daley, 1980), and regardless of the 
weighting, this component has no influence on the rotational 
component . 

As a further demonstration of the variational method, 
the same analysis was balanced two different ways. In one, 
the mass variables were weighed very heavily (W<£> = 1000) 
poleward of +30° latitude. In the other, they were given no 
weight. The 500 mb vorticity and geopotential fields are 
shown in Figs. 9 and 10 with the analyzed contours for these 
two cases. From these figures, it is clear that this range 
of weighting is sufficient to map exactly either vorticity 
or geopotential onto the rotational manifold. This is 
probably not a good method for computing mass from motion 
and vice versa, because the relationships are linear. In 
fact, winds computed from the method compare less favorably 
with observations than the 12-hour forecast first guess. 

To correct this problem, Phillips (1982a) proposes a method 
whereby the nonlinear component is removed from the analysis 
prior to balancing. 

The variational method described in this section is too 
cumbersome to impose full variability into the weights. As 
a result, restrictions that only allow variations with 
latitude were imposed. Such limitations are only necessary 
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Fig. 9. Updated 500 mb vortlcity when (a) weight on mass 
is zero, (b) weight on mass is 1000 poleward of 
30°. The unprocessed analysis is dashed and the 
contour interval is 25 * 1 0 sec -1 . 
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Fig. 10. Updated 500 mb geopotential when (a) weight on 

mass is zero, (b) weight on mass is 1000 poleward 
of 30°. The unprocessed analysis is dashed 
and the contour interval is 60 m. 
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in the normal mode method. As shown in the next section, 
the variational method that uses the balance equation as a 
constraint is capable of using weights with large spatial 
variation. 
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IV. THE BALANCE EQUATION INITIALIZATION 



A. THE VARIATIONAL FORMULATIONS 

For simple models, the nonlinear normal mode balance 
is nearly equivalent to applying the balance and omega 
equations as constraints to the initial conditions (Leith, 
1980). In the variational method developed by Sasaki (1958, 
1970), Stephens (1972) and Haltiner et al. (1976), and 
described in this section, the balance equation is utilized. 
But, instead of the omega equation, the vertical motion is 
derived from the forecast first-guess divergence. 

To impose the balance equation as a constraint, the 
functional 

IU.'M * JaU-^q) 2 + 3 ( W- W Q ) 2 + 2xg m( , <j> ) dA (4.1) 
is minimized. Here, 4> is geopotential, W is the vector 
wind, ^ is stream function, A is the horizontal area over 
which the integral is applied, and Xg is the Lagrange 
multiplier. The constraint is 

m( ^ , 4 > ) = v-f ?ip + 2 J( u , v) - v 2 <j> = 0, (4.2) 

where u and v are the east and north components of wind, 
respectively, J is the Jacobian operator and all operators 
are applied in spherical coordinates. The minimization is 
achieved when the first variation of I(<j>,^) with respect to 
$,Xg and i p is set to zero. Neglecting the variation of 

J(u,v) and noting that 

J v-( ) dA = 0 (4.3) 

A 
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(4.4) 



for integration over the globe gives 

<t> = l + v 2 x b 

and 

SV ip 3 V8*vU-^) + 8V > + vf-VXB + ^V^Xb- (4.5) 

The Eu 1 er-Lagr ange equations are (4.2), (4.4) and (4.5) and 
the unknowns are ip , <|> and Xb* The solution of the Euler- 
Lagrange equations is accomplished by eliminating ^ and <j, in 



using (4.4) and (4.5). This gives: 




? 4 Ux B ) n - f v 2 (AX B ) n = mU"' 1 ,^ n_1 ) , 


(4.6) 


( A<j> ) n = v 2 (AX B ) n , and 


(4.7) 


_2 / . . * n fv 2 (4X B ) n + vf-vUx B ) n 

7 ( A^ ) = 


(4.8) 



The superscript, n, is the iteration count and the delta 
specifies the difference between the current and previous 
i ter at i on , i . e . , 

(A*) 0 = - <D n_1 , (4.9) 

where 

<j> n "^ = <j> when n=l . (4.10) 

Note that all boundary conditions are automatically 
eliminated when the solution is desired for a full sphere. 

The second term on the right of (4.8) was dropped because it 
caused unrealistic adjustments to the winds near the 
equator. Consequently, only elliptic equation solutions for 
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p 

v (ax) and are necessary to solve (4.6)-(4.8). Each 
iteration through (4.6) -(4.8) reduces the residual of the 
balance equation constraint by an order of magnitude, 
therefore two cycles produce sufficiently balanced results. 

A direct solution technique for elliptic equations was used 
to solve (4.6) and (4.8) (see Swarztrauber and Sweet, 1975; 
Rosmond and Faulkner, 1976). 

The forecast first-guess divergence is retained at the 
update time in the following manner: First, the variables 

are balanced as described above. Secondly, the rotational 
component of the forecast first-guess wind field is computed 
and subtracted from the balanced wind. Finally, the 
variables are interpolated to model coordinates where the 
wind, which is now a nondivergent correction, is added to 
the forecast first-guess wind field. The mass variables, 
however, are balanced and interpolated to model coordinates 
in the conventional manner. In this way, the erroneous 
divergence near steep terrain regions remains small, and the 
first-guess divergence is untouched. 

Another method of retaining the forecast first-guess 
divergence is to balance the corrections rather than the 
updated values. This makes the nonlinear term a little more 
cumbersome, but this procedure has the advantage of not 
affecting areas without new data. In this case, the 
nonlinear term is defined as 

J' (u, v) = J(u,v) - J(u 0 ,v 0 ) , 



( 4 . 11 ) 



where the prime is used to designate a correction value 
rather than an updated value. The integral to be minimized 
i s 



where the constraint is 

mU ' ><f> ' ) = v-fvifi' + 2 J ' ( u , v ) - v 2 <j>' = 0. 

This basic approach has also been proposed by Phillips 
(1977). 

B. VERTICAL FILTERING WITH EMPIRICAL ORTHOGONAL FUNCTIONS 
The major problem with the foregoing formulations is 
that the temperature adjustments are not small everywhere. 
For example, if balancing caused the 925 mb height to 
change 10 m and the 1000 mb height did not change, the 
resulting temperature correction between these two levels 
would be 4°C. Such inconsistencies are difficult to avoid 
for grids covering very large regions. 

To minimize the effects of these inconsistent modifica- 
tions in the vertical, the variables are vertically coupled 
by projecting them onto basis functions prior to the 
balance. Because filtering requires that some of the 
unimportant components be dropped, the empirical orthogonal 
functions are the most suitable functions f.or this purpose. 
These functions are derived so that they form a basis set 



I ( <p ' * ^ ' ) 
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that produces the least error when a partial series is used 
to describe a particular three-dimensional field, such as 
geopotential (see Appendix C) . 

An example of the efficiency of the empirical 
orthogonal functions is shown in Fig. 11 for a ten-level 
analysis of geopotential. Here, the first mode explains 95% 
of the total variance, whereas the first two modes explain 
98%. Only four modes are necessary to explain 99.9% of the 
total variance. Projecting the wind analysis onto the first 
four modes retained the details of the jet streams produced 
by the analysis. Holmstrom (1963) first noticed this very 
rapid convergence of the empirical orthogonal functions in 
representing geopotential profiles. The smoothness of the 
first four empirical orthogonal functions insures that 
inconsistent vertical variations of wind or geopotential 
between adjacent levels will be truncated when used to form 
a partial series of the analyzed fields prior to balancing. 
Another advantage of this method is that considerable 
computer time is saved. Rather than treating 10 levels, 
only the four vertical modes need to be balanced. 

To show how empirical orthogonal functions are 
incorporated into the variational balance, the balance 
equation is written in vertical vector form 

m(i|>,$) = v • ( f v tj») + A - v 2 $, (4.14) 
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FIRST 



SECOND 







Fig. 11. The first six empirical orthogonal functions 
derived from a ten-level analysis of 
geopotential. 



70 



where the underlying tilde signifies a column vector, such as 



0 = 



*io( A > 9 ) 



<f>gU,e) 



<J>1 ( * > 9 ) 



(4.15) 



Here a is longitude, e is latitude and the subscript is the 
level number. A is the Jacobian given in (4.2). 
Multiplication of (4.14) by one of the empirical orthogonal 
functions, E-j" 1 ", gives 

m i = vf • Vip 1 + A 1 - v 2 <p 1 = 0, (4.16) 

and the variational integral imposing this constraint is 

-4> 0 i ) 2 + S( W i - \V 01 ) 2 + 2 Am n d A (4.17) 

A 

The solution of (4.17) exactly follows that for (4.1). 

After (4.17) is solved for N of the most significant 
empirical orthogonal functions, the balanced fields are 
constructed from 



N 

On = Z <Pi E. (4.18) 

~ B i=l 1 - 1 

and 



N 

VX \V R = 2 ( v X W • ) E . (4.19) 

— B i=1 i .i 



As shown by (4.19), the balanced vorticity is constructed 
from the empirical orthogonal functions, and then the 
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balanced winds are computed from vorticity. The degree of 
filtering is governed by the size of N, which was set to 
four to keep the jet streams sufficiently strong. 

The filtering described in this section greatly 
improved the erroneous temperature problem. A temperature 
adjustment profile for a filtered and unfiltered balance is 
shown in Fig. 12. As has been typical, the lowest layers 
were particularly poor in the nonfiltered balance, where 
adjustments were in excess of 3°. Zonal averages of the 
temperature adjustments of the lowest layer in the tropics 
were -3°C to -4°C. The filtered balance, on the other hand, 
produced adjustments that were much smaller. The zonal 
averages of the adjustments were everywhere less than 0.5°C, 
and individual adjustments in the lowest layers were about 
1°C. 



C. WEIGHT STRUCTURE 

Typically, the weight a variable receives during the 
variational balancing depends on an intricately computed 
error structure function, which can be derived from optimum 
interpolation methods such as given by Gandin (1963). Such 
error structures are not easily derived from successive 
correction methods, but the amount of data that influences 
any point can be used to Infer analysis accuracy to some 
extent . 
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Fig. 12. Adjustments made to the temperature in order to 
balance wind at 60°E, 30°N on 27 Mar 1982. 

No vertical coupling (dashed line) and vertically 
coupled using empirical orthogonal functions 
( so 1 l d 1 i ne). 
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In the methods described in this section, the weight on 
the geopotential is unity and the weight on the winds is 8, 
thereby keeping the weight function a single variable. 
Limitations to 8 were necessary in the region bounded by 
+30° latitude. Here, g had to be at least 10 4 m 2 sec -2 to 
avoid convergence problems. Outside this region, no such 
'limitations were found to exist. 

The wind weight has two basic parts. One part is 
purely a function of latitude and is prescribed without 
regard to observation density. For this part, the weight 
was set at 40,000 m 2 sec -2 over the region bounded by +30° 
and at 4,000 m 2 sec" 2 over the remainder of the globe. This 
means that a 5 m sec"* modification in wind would correspond 
to a 100 m or a 30 m change in height, depending on where 
the change took place. The second part of the weight 
depends on the number of observations used to describe mass 
or motion, and has a range of values from -2000 m 2 sec -2 to 
2000 m 2 sec" 2 . For example, if a motion variable were 
influenced by five or more observations and there were no 
mass observations, the second part would have a value of 
2000 m 2 sec" 2 . But if the reverse were true, i.e, the mass 
variables were supported by five or more observations, then 
the second part would have a value of -2000 m 2 sec" 2 . The 
final weight value is the sum of the first and second parts. 
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In this chapter, the classical balance equation has 
been incorporated into a global initialization that uses 
calculus of variations. The method involves the solution of 
elliptic equations on a sphere, but otherwise the method is 
simple and flexible in terms of variable weighting. This 
system was thoroughly tested by Barker (1981a). In the next 
section, the role of the balance equation in data assimila- 
tion will be compared to that of the nonlinear normal mode 
method . 
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V. COMPUTATIONAL RESULTS 



A. EXPERIMENT DESCRIPTION 

The results comparing the different initialization 
methods were obtained using the FGGE level 1 1 a data set for 
the period between 16 Nov 1979 and 21 Nov 1979. The level 
1 1 a data were available on an operational basis, conse- 
quently no special processing was performed for this set. 

The observing systems available at this time are described 
by Fleming et al. (1979). 

Four different initialization methods were tested in 
data assimilation tests that covered the period of the data 
set. Two of the methods tested used some form of the 
balance equation, and the other two used the variational 
nonlinear normal mode method. 

In the two balance equation methods, the divergence was 
obtained from the forecast first guess. In the first 
method, hereafter referred to as BE1, the corrections to the 
forecast were balanced and interpolated to the model coordi- 
nate. This maintained a balance between mass and motion and 
preserved the divergence that was present in the forecast 
just prior to the time of the update. The second method, 
which shall be called BE2, balanced the total analysis 
(forecast first guess plus corrections) and then converted 
the balanced winds to nondivergent corrections. The motion 
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variables were then treated as in B E 1 , except the mass 
variables were treated as the new values in the model. The 
main difference between the two methods is that the correc- 
tions were balanced in BE1 and the updated variables were 
balanced in BE2. Both methods used the same weight assign- 
ment described in Chapter IV. 

Two variations of the normal mode approach were used. 
In the first case (NM1), the weights were varied with 
latitude. Poleward of +30 degrees, the weight on the 
geopotential was 10, that is an order of magnitude larger 
than naturally occurs in the normal mode method. Further 
increases in geopotential weight make only slight improve- 
ments in the fit of geopotential, as discussed in Chapter 
IV. 0. Equatorward of +30 degrees, the geopotential weight 
was 0.5. In the second case (NM2), the geopotential weight 
was unity everywhere. 

B. DIFFERENCES BETWEEN BALANCED AND ANALYZED VARIABLES 

The objective analysis method described above is used 
to fit the available observations relative to the first- 
guess fields from the forecast model. Unfortunately, large 
regions of the earth have inadequate data coverage in terms 
of quantity and quality. Mass corrections are frequently 
made in regions without motion corrections and vice versa. 
This, in turn, causes large imbalances between the mass and 
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wind fields which often cause the model to reject much of 
the analyzed information (Daley, 1980; Daley and Puri, 

1980). Similarly, the initialization may reject updated 
information by balancing to the motion variables when the 
mass variables are more correct. To show how the balanced 
conditions differ from the analyzed ones, the differences 
between the balanced and analyzed variables were plotted for 
the different balancing methods. 

The RMS averages of the analyzed corrections and the 
differences between the balanced and unbalanced conditions 
are shown in Fig. 13 for 8E1 and 8E2. Except for tempera- 
ture, the two methods produced about the same modifications 
to the analyses. Surprisingly, these changes to the 
analyses were nearly as large as the original corrections. 

The globally averaged temperature modifications are 
much larger for the BE2 method than the BE1 method, 
primarily because of the adjustments made to the lower 
levels of the model. These adjustments were much smaller in 
subsequent update cycles ( 0.8 degrees), but were consis- 
tently twice as large as the BE1 method. 

The large differences between analyzed and balanced 
winds are related to the inability of the objective analysis 
to project the corrections onto the rotational component of 
the wind. This is further illustrated in Fig. 14. 
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Fig. 13. RMS differences between the first guess and analysis (solid line), 
and RMS differences between balanced values and analysis (dashed 
lines) where BE1 is represented by long dashes and BE2 is represented 
by short dashes. The subscript 4 designates level 4, and the global 
RMS values are given above the respective plots. 
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Fig. 14. Divergence (a) and vorticity (b) at 500 mb in the region of New 

Zealand. The values after balance are depicted with a solid line and 
the a n a 1 y 7 e d values are depicted with a dashed line. Contour interval 
is 10 10"6 sec -1 for divergence and 25 10"6 sec -1 for vorticity. The 
time is 00Z GMT 11 Nov 1979. 



Regions such as the one shown have few wind observations. 
When there is an isolated report, it produces a correction 
that is largely divergent (see Barker, 1981b), and regard- 
less of the weighting on the winds, this divergence is 
rejected during the balance. The rotational component, 
however, can be fit as accurately as desired (see Chapter 
IV. D), by adjusting the wind weights. The vorticity 
shown in Fig. 14 is reasonably close to the analyzed 
values considering the moderate values of weighting applied 
to the winds in this area, (j^-)^ ' 4000 m^ sec"^. 

The RMS average differences between normal mode 
balanced and unbalanced variables were also compared to the 
RMS average of the analyzed corrections (see Fig. 15). In 
these comparisons, however, the divergent component of the 
wind correction was removed from the analysis before inter- 
polation to the model coordinates. The method with stronger 
weight on geopotential produced closer fits of temperature 
and pressure than did the other method, but it failed to fit 
the winds as well. This is to be expected in the varia- 
tional method. The globally averaged wind differences were 
between 1.0 and 1.5 m sec"-*-, but when divergence was not 
removed from the analysis, the differences were consistently 
larger (between 2.0 and 2.7 m sec"^). The surface pressure 
differences were as large as the analyzed corrections over 
much of the earth. 
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Comparing the balance equation results (Fig. 13) to 
those of normal mode methods (Fig. 15) indicates that both 
techniques were similar in the way that they modified 
analyzed variables. The balanced winds differ from the 
analyzed winds by about the same amount for all methods when 
the analyzed divergence was not removed prior to balancing 
(these curves are not shown). The normal mode methods were 
slightly better at fitting the analyzed temperature, but 
poorer at fitting the surface pressure. 

The modifications to the analyzed variables were 
significantly large in all experiments performed. However, 
whether the impact of the balancing can be considered 
beneficial to the data assimilation process principally 
depends on the magnitude of the gravity wave noise that 
still exists. In the next section, the noise levels that 
are present in the model before and after balancing are 
compared for the balance equation and normal mode methods. 

C. ELIMINATION OF GRAVITY WAVE NOISE 

In this section, gravity noise produced in forecasts 
initialized with no balance, the balance equation and the 
normal mode method are compared. In these comparisons, 
surface pressure tendency was used as a measure of the 
gravity noise. Although this quantity reflects only the 
presence of external gravity waves, the highest frequencies 
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clashes) and NM2 (short dashes). 



come from these external waves. Except where stated 
otherwise, a dry version of the global model was used to 
perform the integrations. 

In the first comparison, the effectiveness of the 
balance equation was tested by comparing it with an 
initialization that just removed divergent winds from the 
update corrections. In the balance equation method used for 
testing, the corrections were balanced, interpolated to 
model coordinates and then added to the forecast first guess 
(method BE2). As can be seen in Fig. 16, the balance equa- 
tion forecast is less noisy than the forecast initialized 
with nondivergent winds, particularly during the final hours 
of the forecast. The slower adjustment rate of the 
unbalanced forecast is probably due to the scale of the 
gravity waves, since global models may carry large gravity 
waves that do not readily disperse their energy (Bourke, 
1972). This version of the balance equation initialization, 
on the other hand, is effective at removing these large 
scale waves, but it does not balance around terrain. 
Consequently, the balance equation method contained notice- 
able small scale noise. 

Using the model's normal modes, the linear balance was 
compared to the nonlinear balance (Fig. 17). The nonlinear 
balance forecast required one hour adjustment time and 
resulted in about one half as much noise as did the linear 
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Fig. 16. Global RMS average of surface pressure tendency from a forecast 
initialized with no divergence (solid line) and a forecast 
initialized with the balance equation (dashed line). 
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Fiu. 17. Global RMS average of surface pressure tendency from a forecast 

initialized with linear normal mode balance (solid line) and nonlinear 
normal mode balance (dashed line). 



balance forecast. These results illustrate the effective- 
ness of the nonlinear balance step, which can balance around 
terrain (Daley, 1979; Williamson and Temperton, 1981) and 
can include the nonlinear component of the balance 
( Machenhauer , 1977). 

Plots of surface pressure tendencies from predictions 
initialized with the balance equation and with the nonlinear 
normal mode method are shown in Fig. 18. These results show 
the superiority of the normal mode method over the balance 
equation. After seven hours, however, the forecast started 
from the balance equation contained only slightly more noise 
than the forecast from the normal mode method . 

Finally, the impact of adding latent heating effects to 
the forecasts is illustrated in Figs. 19 through 21. The 
increased noise level is most noticeable in the normal mode 
runs (Fig. 20) where the forecast seems to have required 
about six hours of adjustment time. The balance equation 
run was only slightly affected by the heating (Fig. 19), 
possibly because less precipitation was generated during the 
early hours of the forecast. In any case, the result is 
that the normal mode method was no better at noise 
suppression than the balance equation when latent heating 
effects were included in the forecast . This result was 
observed consistently several times during the course of 
development of the methods. 
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Fig. 18. RMS of surface pressure tendency for the normal mode initialization 
(dashed line) and for the balance equation (solid line) when latent 
heating is not included in the forecast model. 
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Fig. 19. RMS average of surface pressure tendency in forecasts initialized with 
the balance equation method. The adiabatic forecast is represented by 
a dashed line and the forecast including latent heat effects is 
represented by a solid line. 
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Similar to Fig. 19 except for the nonlinear normal mode technique. 
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Fig. 21. RMS of surface pressure tendency for the nonlinear normal mode 

Initialization (dashed line) and for the balance equation (solid line) 
when latent heating is Included in the forecast. 



In summary, the normal mode method was most affected by 
latent heating in the forecast. This was probably because 
it allowed more precipitation to occur during the early 
hours of the forecast compared to the balance equation 
method. Considering that one problem associated with 
initialization has been the lack of precipitation early in 
the forecast, the noise level increase in the forecast 
including latent heating may be a symptom of a beneficial 
result. It is noteworthy that latent heating effects 
generated twice as much gravity noise as the integrations 
that did not include latent heating. This was true even 
after the initial adjustment period was complete. To 
explore the question of latent heating effects more fully, 
comparisons of the precipitation, vertical motion and 
cyclogenesis during data assimilation runs were made. The 
results from these runs are described in the following 
section . 

0. VERTICAL MOTION PRECIPITATION AND CYCL06ENESIS 

The lack of precipitation early in the forecast is 
considered to be a major problem in the initialization of 
numerical models (Tarbell el al., 1981; Bengsston, 1981; 
Daley, 1981). This is particularly true of forecasts 
initialized without vertical motion. An example of the 
problem is shown in Fig. 22. These two forecasts were 
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initialized with the balance equation with no provisions to 
include vertical motion. The forecasts for the region 
west-northwest of the Hawaiian Islands predicted large 
precipitation rates after six hours, but only slight amounts 
before this time. 

These initially small precipitation rates occur when no 
vertical motion is included in the initialized data. Inclu- 
sion of the omega equation (Tarbell et al., 1981), nonlinear 
normal mode initialization (Leith, 1980) and retention of 
the forecast first-guess vertical motion provide possible 
solutions to this problem. Unfortunately, the omega equa- 
tion and normal-mode methods have not worked well in the 
tropics (Bengsston, 1981; Tribbia, 1981), and the forecast 
first guess may suffer from inaccuracies in the forecast. 

The balance equation methods discussed in Chapter III 
use the forecast first guess divergence. The normal mode 
initialization partially recomputes this divergence through 
the nonlinear balancing of the external and first internal 
vertical modes. The forecast and derived divergence from 
the normal mode methods are given in Fig. 23. A sequence of 
forecast and computed divergence fields is given in Fig. 24. 
Notice tnat only small differences exist between the fore- 
cast and computed divergence patterns. This similarity 
between the two divergent winds suggests that the forecast 
divergence is a fairly accurate quantity. 
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Fig. 23. The first-guess divergence (dashed line) and 
normal mode computed divergence (solid lines) 
at 500 mb valid at 00 GMT 20 Nov 1979. The 
contour interval is 10" 5 sec" 1 . 
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Fig. 24. Forecast first-guess divergence (a, c, e) and 

divergence computed using the normal mode method, 
NM2 (b, d, f), for three successive 12-hourly 
updates beginning at 12 GMT 17 Nov 1979. Contour 
interval is 10" 5 sec" 1 . The interval between 
- 1 - 1 0 “ 5 sec" 1 and -2*10" 5 sec" 1 is cross-hatched 
and the interval between 10" 5 sec" 1 and 2*10 
sec" 1 is cross-hatched twice. 
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Precipitation rates from five forecasts that are 
identical except for the initialization procedure are shown 
in Fig. 25. In the first case (Fig. 25a), no initialization 
was used. Even the excessive divergence produced by the 
analysis was included in the initial data. The other 
methods, which are designated BE1, BE2, NM1 and NM2, are 
described in Chapter V.A. These forecasts all produced a 
nearly identical precipitation pattern along the coast of 
North America. The normal mode methods (NM1, NM2), however, 
produced less precipitation in the central Pacific than did 
the balance equation methods (BE1, BE2). However, on the 
other hand, normal mode methods have produced the largest 
precipitation rates south of Japan. None of the methods 
produced a persistent bias in the strength of precipitation 
rates. Considering the sensitivity of precipitation rate to 
small changes in initial data, however, the similarity of 
the patterns is quite remarkable. In particular, the lack 
of spurious precipitation in the forecast without any 
balancing is most surprising. 

The precipitation rates during the first twelve hours 
of forecasts initialized by the balance equation and normal 
mode methods are shown in Fig. 26 for a region just north of 
South America. Separate rates are given for the first and 
second six-hour periods to illustrate the impact of initial- 
ization. During the first six hours of the forecast, the 
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Fig. 25. Precipitation rates during the first six hours from forecasts 

initiated from: (a) no balance, (b) BE1, (c) BE2, (d) NM1 and (e) NM2. 
The contour interval is 2 cm h r and the interval between the 6 and 8 
isolines is cross-hatched. The time is 00 GMT 16 Nov 1979. 
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Fig. 25 . Continued . 
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Fig. 25. Continued. 




Fig. 26. Precipitation rate during the first six hours from 
forecasts initialized with BE1 (a), BE2 (c), NM1 



(e) and NM2 (g), and precipitation rate during the 
second six hours for BE1 (b), BE2 (d), NM1 (f) and 
NM2 (h). The contour interval and cross-hatching 
are as in Fig. 25. The starting time is 00 GMT 11 
Nov 1979. 
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precipitation extended over much larger regions and was 
stronger than during the second six hours. This tendency is 
just the opposite effect observed in the nondivergent case, 
and helps to explain why the model produces more gravity 
wave noise when the latent heating effects are included in 
the forecast. The precipitation patterns appear to be less 
similar between the forecasts during the second six hours 
than during the first six. In fact, after several update 
cycles, the precipitation fields produced by the different 
methods contain little similarity in the tropics. 

Slight differences due to the various initialization 
methods may cause the assimilation runs to diverge slowly 
from each other. The largest differences are likely to 
occur in regions of strong baroclinic instability where 
precipitation can play a role. To examine ?his possibility, 
a moderately intense surface low, which developed along the 
Aleutian Islands, was studied. This case of surface cyclo- 
genesis of 20 mb in 24 hours was supported by an upper level 
short wave that traveled along the island chain in twenty- 
four hours. Three twelve-hour forecasts are shown for the 
various initialization methods tested (Figs. 27 through 34). 
Comparing forecasts rather than analyses helps to eliminate 
the possibility that gravity modes allow a closer fit to the 
data than actually exists by the meteorological modes. As 
can be seen from the verifying analyses (Fig. 35), none of 
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Fig. 27. Twelve hour forecasts of 500 mb geopotential (a, 
c, e) and sea-level pressure (b, d, f) during the 
assimilation run using the 3E1 initialization 
method. Contour interval for geopotential is 60 m 
and for sea-level pressure is 4 mb. The 4920 to 
4980 m interval is cross-hatched in the 500 mb 
maps. The 996 to 1000 mb and 980-984 mb 
intervals are cross-hatched in the sea-level 
pressure maps. 
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Fig. 28. Precipitation rate during the first six hours (a, 
c, e) and second six hours (b, d, f) of the 
forecasts from the BE1 initialization. Contour 
interval is 2 cm hr -1 and the contour 
interval between 2 and 4 cm hr" 1 is 
cross-hatched. 
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Fig. 30. Similar to Fig. 28 except for the BE2 method. 
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Flg * 31. Similar to Fig. 27 except for the NM1 metnod. 
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Fig. 34. Similar to Fig. 28 except for the NM2 method 
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Sea-level pressure analysis for 19 Nov 00 GMT produced by the BE1 (a), 
BE2 (b), NM1 (c) and NM2 (d) assimilation runs. 



the forecasts developed the system as rapidly as the 
atmosphere. However, the normal mode method came closest by 
developing slightly deeper systems each time. The NM2 
method managed to do slightly better than the NM1 method. 

The weakest system was produced by the BE2 method, but it 
was only 4 mb too weak. The normal mode methods generated 
the largest precipitation amounts. This is particularly 
true during the 6 to 12 hour forecast made from the 17 
November 1979 12 GMT data. During most of the period, 
however, the precipitation differences were quite small. 

In general, forecast experiments that used divergence 
in the initial conditions produced precipitation rates in 
the early hours of the forecast that were larger than those 
produced after six hours. This is the opposite effect 
observed in forecasts using a nondivergent initialization. 
Whether the divergence was computed from the normal mode 
method or derived from the forecast first-guess seemed to 
make little difference in the resulting precipitation and 
divergence. During data assimilation experiments lasting 
several days, however, these differences between the methods 
became more pronounced. A study of cyclogenesis over the 
North Pacific showed that while the normal mode methods 
tended to deepen the system faster and somewhat more 
accurately, the balance equation methods were nearly as 
ef f ect i ve . 
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Unfortunately, these results which describe a single 
case are hardly conclusive. To extend the verification of 
the different initialization methods over broader areas and 
more cases, comparisons of many forecasts to observations 
are made in the next section. 

E. FORECAST VERIFICATION 

Justification for balancing during data assimilation 
runs is readily demonstrated with the verification compari- 
son of a balanced and unbalanced forecast. Forecast 
verifications against all observations for these runs are 
plotted versus latitude in Fig. 36 for geopotential and 
winds. The results indicate that although the wind 
verification is unaffected by the presence of gravity wave 
noise, the geopotential verification is drastically 
affected. This forecast from an unbalanced initial state 
has RMS errors almost double that of the balanced forecast 
over much of the earth. Such large errors in a forecast 
first guess may cause the quality control programs to reject 
observations that should not be rejected and to add noise to 
the resulting analysis. (Heavy damping filters applied 
during integration may make the forecasts more usable, 
however . ) 



>500 




LATITUDE 



Fig. 36. Latitudinal variation of RMS (a) pressure height 
and (b) wind differences between observations and 
12-hr forecasts without balance (stars) and with 
NM1 initialization (circles). 
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While it is certain that the forecasts from an 
unbalanced initial state are unsuitable as a first guess for 
an analysis of geopotential observations, and therefore are 
not useful in data assimilation, the most effective 
balancing procedure to be used is not clear. Some small 
differences do exist between the assimilation results using 
different initialization methods. This can be seen (Fig. 

37) in the analysis of a short wave upstream of Australia. 
The 12-hour forecasts (rather than analyses) are shown for 
this system so that data resolved onto gravity waves could 
be dispersed by the model. The normal mode methods (NM1 and 
NM2) produced a slightly deeper wave than did the balance 
equation methods (BE1 and BE2). This is consistent with the 
results in Chapter V.C. 

An extensive objective verification study was performed 
for the different methods, where the RMS differences between 
the 12-hour forecasts and observations were computed. These 
computations were made for ten forecasts covering the period 
of the data assimilation experiments. Unfortunately, this 
type of test may tend to favor the smoother forecasts, and 
therefore the interpretation of results should be made with 
this in mind. 

Although verification was performed by region as well 
as globally, the regional verification added no new informa- 
tion not readily apparent in the global statistics. The 
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Fig. 3 7. Sequence of 12-hr forecasts of 500 mb height from assimilation runs 
using BE1 (a), BE2 (b), NM1 (c) and NM2 (d) initialization methods. 
The contour intervals 4860-4920 m and 5760-5820 m are cross-hatched. 



verification is presented by observation type, which tends 
to regionalize the verification to some extent. For 
example, observations from ships, satellites and aircraft 
are primarily taken over the oceans, whereas radiosonde 
observations are mostly taken over land. 

Forecast verification against radiosonde geopotential 
and wind observations are presented in Figs. 38 and 39 for 
the four methods of initialization. The lowest RMS differ- 
ences between forecasts and observations occurred in the 
balance equation method, BE2. However, the differences 
between methods are very small. For geopotential, the 
differences produced by the various methods are generally 
less than 2m, which is only about 4% of the total RMS error. 
For wind, this difference is generally less than 3 m sec - ^, 
which also represents about 4% of the total error. 

Comparing the normal mode methods, it can be seen that NM1 
verified against geopotential observations slightly better, 
whereas NM2 verified against wind observations slightly 
better. This reflects the emphasis that the variational 
balancing placed on the respective variables. However, 
notice that 8E1 verified geopotential as well as the NM1 
method and wind as well as the NM2 method. 

Forecast verifications against ship observations of sea 
level pressure, which were converted to geopotential through 
tne hydrostatic equation, are given in Fig. 40. The 
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rig. 38. RMS differences between radiosonde observations 
of pressure height and 12-hr forecasts for the 
assimilation runs comparing B E 2 (a), NM1 (b) and 
NM2 (c) with BE1 initialization methods. Each 
data point represents the error in the assimila- 
tion model just prior to the update. Abscissa 
labels are hours after start of the assimilation 
run . 
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c ig. 39. Similar to Fig. 38 except for radiosonde wind 
observations. 
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Fig. 40. Similar to Fig. 38 except for surface ship 

observations converted to geopotential height. 



improvement of the BE2 method over the other methods is much 
more significant for this data. Compared to BE1, the 
improvement is as high as 10% of the total RMS error. The 
NM1 method, which weighted the geopotential analysis more 
than the NM2 method, verified somewhere between BE2 and BE1, 
whereas NM2 produced approximately the same results as BE1. 

Verification using satellite-generated geopotent 1 al s 
(Fig. 41) again show the BE2 method to be superior to the 
other methods. The improvement ranges between 6 and 12%. 

The other initialization methods resulted in rather similar 
verifications. 

Satellite wind forecast verification (Fig. 42) shows 
that the BE2 method gave about 0.5 m sec'^ smaller forecast 
error, or about a 5% improvement over the other methods. 

This verification is based mainly on the low level ( 925 mb) 
satellite observations in the tropics. Once again, the 
other verifications are similar to each other. 

Unlike the other types of data, aircraft wind forecast 
verification failed to show much difference between the 
methods (see Fig. 43). Since aircraft observations are 
primarily taken between 300 and 250 mb, this suggests that 
the four methods produced comparable results at the upper 
levels over the oceans. 
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Fig. 41. Similar to Fig. 38 except for satellite derived 
geopotent i al s 
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Fig. 42. Similar to Fig. 38 except for satellite 
observations . 
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Fig. 43. Similar to Fig. 38 except for aircraft 
observat ions 
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wind 



In summary, the BE2 method, which uses the balance 
equation to var i ati onal ly balance the complete analysis 
(forecast first guess plus corrections), produced the best 
verification scores. This result is most pronounced at the 
low levels over the oceans. Visual inspection of the fore- 
casts produced by this method indicate that they were also 
smoother than forecasts produced by other methods. While 
RMS scores tend to be better for smooth fields, the much 
better verifications against ship observations produced by 
the BE2 method are difficult to dispute. The differences 
among the remaining methods are too small to conclude which 
is best. 
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VI. SUMMARY AND CONCLUSIONS 



In this study, two different static initialization 
methods have been developed and tested in data assimilation 
experiments. These methods not only reduce gravity wave 
noise, but they also fit the meteorological modes to the 
most accurately analyzed variables through calculus of 
variations. The methods are based on the balance equation 
and a nonlinear normal mode procedure. Of the two methods, 
the normal mode method would be the more cumbersome to 
apply, because the most general form with the weighting 
fully variable in the horizontal requires more computer 
memory than was available. 

The two methods were compared in various ways. The 
comparison tests were designed to show how the analyzed 
fields were modified during the balancing, how each 
controlled gravity wave noise, how precipitation varied 
during the early hours of the forecast, how the forecast of 
a rapidly developing system was affected, and finally, how 
short range forecasts from the various methods verified 
against observations. 

Both methods modified the analyzed fields during the 
balance by a large amount when compared to the size of the 
corrections (differences between analysis and forecast first 
guess). Surface pressure modifications were particularly 



large, sometimes exceeding the magnitude of the corrections. 
Modifications to the winds were primarily caused by 
unrealistic divergence patterns in the analyzed winds, which 
were nearly as large as the analyzed corrections. When the 
divergence was removed prior to balancing, the modifications 
were much smaller. For wind and surface pressure, the 
magnitude of the adjustments was very nearly the same for 
the two methods. However, the modifications to the tempera- 
ture analysis were larger for the balance equation method 
when it was used to balance the total analysis (corrections 
plus first guess) than when it was used to balance just the 
corrections . 

In terms of gravity wave noise removal, the normal mode 
methods performed significantly better in a dry version of 
the prediction model than did the balance equation method. 
Adding the effects of latent heating, however, tended to 
mask this improvement, because the heating effects caused 
gravity wave noise to more than double relative to the dry 
version of the mode 1 . 

The increased gravity noise due to the heating may be 
caused by many factors. It is partly caused by the way in 
which the heating effects are included. Each fifth time 
step, the heating effects are added. A Matsuno time step is 
then used prior to resuming the leapfrog time stepping. The 
gravity wave amplitudes produced are larger than if the 



127 



heating effects are added incrementally over each time step. 
Another factor that may make the latent heating effects more 
pronounced is that the model tends to be warmer and drier in 
the lower layers than observed in the atmosphere (Johnson, 
1976; Payne, 1981). After each update, the solution tends 
toward the model equilibrium state. Also, because no 
analysis is made of moisture, changes in the mass and motion 
structure may require that the model undergo some adjustment 
before latent heating matches the updated systems. None of 
these effects are controllable by the initialization proce- 
dures under study, however. 

Integrations performed without divergence (vertical 
motion) required several hours to develop realistic 
precipitation rates. This was particularly true in the 
tropical regions in a test with the balance equation 
initialization that did not include divergence from the 
forecast first-guess. 

Comparisons of divergence from the forecast first-guess 
and that computed using the nonlinear normal mode balance 
applied to the external and first internal vertical modes 
revealed only small differences. Other comparisons of 
precipitation forecasts from the balance equation, normal 
mode and no balance conditions showed only minor differ- 
ences. In the balance equation and no balance methods, the 
forecast first-guess divergence was present in the initial 
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conditions. The unbalanced forecast surprisingly produced 
no noticeable spurious precipitation forecasts, even when 
the unrealistic divergence produced during the wind analysis 
was included. Unlike the forecast initialized without 
divergence, however, largest precipitation rates were pre- 
dicted during the first few hours of the forecast. It 
appears, then, that whether the forecast first-guess 
divergence is partially recomputed using the normal mode 
methods or even mixed in with unrealistic divergence, little 
difference will exist in the precipitation forecast. 

The effectiveness of the model in assimilating data 
around a rapidly developing surface low over the Pacific was 
examined for the different initialization methods. The 
representation of the cyclone development was very similar 
for the various methods. The maximum central pressure 
difference between the methods was 4 mb. The more intense, 
and slightly more accurate representation of the low was 
produced by a normal mode method, whereas the least accurate 
was produced by the balance equation method that balanced 
the total analysis and used the forecast first-guess diver- 
gence. The results from this single example are only 
suggestive, however. 

To produce results over a wider number of cases, the 
forecast first-guess fields were verified for the four data 
assimilation methods. Two methods used variations of the 



balance equation method and two others used variations of 
the normal mode method. The results from these runs showed 
that while minor differences existed between the forecasts 
for much of the data, the balance equation method that 
balanced the total analysis, rather than the corrections, 
produced the most accurate short-range forecasts near the 
surface over the oceans. The forecasts produced from the 
normal mode methods tended to contain slightly stronger 
systems than the balance equation methods, which may have 
actually reduced all the verification scores for this method 
relative to the smoother fields from the balance equation 
method . 

In conclusion, the results of this study show that it 
is possible to initialize a model with the forecast first- 
guess divergence. This allows continuity in the precipita- 
tion rates produced by the model during the updates. 
Consequently the variational balance equation method 
produced results that are competitive with, and in some 
respects, better than the more elegant nonlinear normal mode 
method. This conclusion is based on precipitation rates 
during the early hours of the forecast, gravity wave noise 
and short-term forecast results. In terms of variational 
weight assignment, the balance equation methods are less 
cumbersome, and thereby allow more flexibility. 



As a note of caution, however, it should be mentioned 
that assimilation systems, such as the ones tested, are too 
complicated to guarantee absolutely that the tests were 
without flaws. For example, the balance equation method 
that scored the highest forecast verification scores was 
also the only method that used a slight variation of the 
interpolation to model coordinates. It interpolated 
updated mass variables instead of corrections. Interpola- 
tion of corrections did not insure that the sea level 
pressure under terrain matched observations corrected to sea 
level. Consequently, regions such as the Himalayas 
contained sea-level features that were not present at 
terrain level. Another factor is the type of assimilation- 
prediction that was used. The version of the UCLA model 
used in the assimilation produces much larger surface 
pressure tendencies than does a dry version of the same 
model. This factor, which masked the benefits of the normal 
mode initialization, may not be so prevalent in future 
generation models. 



APPENDIX A 



DATA ANALYSIS WITH SIMULTANEOUS FILTERING 



The inherent smoothing of a successive correction 
scheme can be determined by neglecting discrete spacing of 
the observations and assuming isotropy in the weight 
specification (Barnes, 1973). Under these conditions, 

F(r) =J " u ( e ) f ( r-e ) d e . (Al) 

-oo 

The convolution theorem allows us to take the Fourier 
transform of this equation or 

F(K) = 5(K)f(K) (A2) 

where K is the wave number equal to 2 tt / x, and X is wave- 
length. The hat is used to show that the variable is 
transformed into wave number space, e.g., 

«(K) =/.„ »(e)e iKe d e (A3) 

Recall that F is a complex quantity whose magnitude is 
represented by 

| F | » VFF* ( A4 ) 

or 

I F | = | w i if | , (A5) 

which means that the spectral response of the analysis may 
be identified by ju(K)j. 

If the analysis is repeated to converge on the data 
more closely, the resulting product will be 



(A6) 



Fy(r) = F x (r) + JT" o , 2 ( e ) [f ( r- e ) -F x ( r ) ]de 
where F]_ is the result of applying (Al) with weight function 



<!)]_. This equation transforms to 

Fj(K) = F X (K) + tS 2 ( K)f(K) - u 2 (K)F 1 (K), (A7) 

or 

Fj(K) - [u! ( K ) + £ 2 (K) - u 2 (K)£ 1 (K)]f (K) (A8) 

and the spectral response of two passes is identified by 

1 ,j0 yI = I^i ^2 ~ ^2 aj ll ( A 9 ) 

Specifying the weight function as Gaussian, where 

oi 1 ( e ) = a-j exp(-e 2 /B 2 ), (A10) 

012 (e) = exp(-e 2 /yB 2 ) , (All) 

simplifies the computation of spectral response, y and B c 



are arbitrary constants and a^ and a 2 are normalizing 
coef f icients, e.g., 

a l = [To, ex P( - e 2 /B 2 ) d £ ]-l . (A12) 

The Fourier transform of (A10) and (All) are 

(K) = exp ( -B 2 K 2 /4 ) and (A13) 

u 2 (K) = exp ( -yB 2 K 2 /4 ) . (A14) 

It is now possible to design the spectral response, wy, by 
appropriate choices for Y and B. 
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The objective analysis for discrete spacing of the 
observations is 

N 

F'(r) • z 0J l(-Sk) f ' ( r " £ k) * ( A 1 5 ) 

k = l 

In this equation and those that follow in this section, the 
upper case letters, such as F', designate gridded fields of 
the variable being analyzed, whereas the lower case letters 
are related to observations that are irregularly located in 
space and time. The vector, specifies the separation in 

space and time between the observation, f^, and the grid at 
point r. The quantity f' represents the difference between 
the observation and the forecast that is to be updated and 
F' is the analyzed correction field that results from A15. 

The weighting function, determines how the observa- 
tions are weight-averaged at each point on the grid. There 
is no upper limit to the size of e, but a practical limit 
occurs when wj is sufficiently small as to make correction 
meaningless. The volume over which computation is performed 
is referred to as the scan volume, and it represents the 
region from which observations are weight averaged to 
produce a value for point r. Equation A15 was first 
developed by 3ergthorsson and Doos (1955). Cressman (1959) 
overcame some of the inherent oversmoothing of this 
technique by incorporating multiple analysis scans, with 



each scan using a progress i ve ly smaller radius. In this 
successive correction method, the second pass uses a weight, 
oj 2> to reduce the remaining discrepancy between the analysis 
after the first pass and the observation, f", i.e., 

N 

F"(r) = E “ 2 ( £ l<)f"(r -ei.). ( A 1 6 ) 

k = l ~ 

The resulting correction is equal to the sum of F' and F". 

Although the Cressman method is usually designed to 
perform three or more passes through the data, Barnes (1964, 
1973) and Stephens (1967) have shown that careful design of 
the weight function may provide adequate results after only 
two passes. The weight functions defi ned by Barnes are 

a* i = a]_ exp(-e 2 /3 2 ) ( AT 7 ) 

and 

uj 2 = ^2 exp(-e 2 / Y B 2 ) . (A18) 

As already shown, an analysis that uses ( A 1 8 ) and 
(A19) to compute the weights, and has a sufficiently dense 
observational coverage, will produce a spectral response 
<2 j(K) = exp(-B 2 K 2 /4) + exp ( -y B 2 K 2 /4 ) 

exp [- ( 1 +y) B 2 K 2 /4], ( A 1 9 ) 

where 3 and Y are arbitrary parameters, and K is the wave 
number ( 2 ^ / x » where a is wavelength). The symbol (~) 
signifies that the variable is transformed into wave-number 
space. It is possible to use this equation to design 



variable filtering characteristics that are dependent on 
observation type in the analysis. The parameter B is used 
to limit the volume of influence of an observation, and the 
parameter Y is used to specify the degree of inherent 
f i 1 ter i ng . 

To design the desired filtering into the analysis, it 
is convenient to define the weight function as the product 
of three functions, i.e., 

co = co H (e H )co v (e v )co t (e t ) (A20) 

where the subscripts H, v and t designate the functions 
describing the horizontal, vertical and temporal dependence, 
respectively. In this form, it is possible to make the 
filtering different for each dimension in the analysis. 

A vertical weighting function that allows ample 

% 

vertical variability while maintaining some vertical 
coupling is desirable- The weight function, co v , in (A20) 
is proportional to exp( -e v /B v 2 ) , where e v = ln(P|</P r ) 
represents the pressure separation between the observation 
and analysis level. The constant 3 V is 0.6, which produces 
a vertical scan radius that corresponds to the positive 
values of the prediction error covariances computed by 
Hallett (see Rutherford, 1976). The vertical filtering 
parameter, y v , is allowed to vary with observation type. 

For satellite soundings it is 0.3, whereas for other data 
types it is 0.3. Thus the satellite sounding corrections 



are very smooth with height, thereby avoiding the problem of 
removing inversions that are formed by the model or are 
present in the other observations (Tracton et al 1980). 

For example, updates from satellite soundings containing 
vertical wavelengths of 0.5 are damped to 20% of their 
original values, whereas this same wavelength in a radio- 
sonde is damped only 50%. 

As is the case for the vertical function, the 
horizontal weighting function is designed to limit the 
influence of an observation to an area that roughly 
corresponds to the positive values of correlation structure 
function. Using the curves produced by Buell (1972), B H is 
then 3.24 grid intervals on the 2.5° mesh used by the 
analysis. To account for the spherical geometry, the hori- 
zontal separation between observation and grid point is 
defined by 

e H 2 ■ (ux ) z ♦ y 2 , (A21 ) 

where x and y are the zonal and meridional distances in grid 
intervals, u is the map factor 

y = max [cos 0 , 0. 5 ], (A22 ) 

and 9 is latitude. The lower limit on y distorts the region 
of influence an observation may have poleward of 60°, but it 
prevents obvious computational problems near the poles. 



Inman (1970) has shown that successive correction 
methods have a tendency to diffuse narrow jets such as occur 
in the upper troposphere. An e 1 1 i pso i da lly shaped weighting 
function with the major axis aligned along the wind 
direction tends to avoid the difficulty. Therefore, the 
weighting function for the wind is computed from 



u> H (e H ) = exp ( m y / B H ^ ) , (A23) 

where 

m y = [0.7 + 3.0 sin(9 v -9 £ ) ] , (A24) 

9 v = arctan (v/u) and (A25) 

9 = arctan (y/x) . (A26) 

£ 



Inclusion of the map factor in (A26) made no impact in the 
analysis and therefore was omitted to save considerable 
computer time. 

. The present assimilation system updates the forecast 
every twelve hours, which means that during any single 
analysis time, there are likely to be data whose observation 
time differs from that of the analysis by six hours. As is 
the case in spatial dimensions, a poorly defined time weight 
function will result in damping and shifting of the small 
scale waves. If the weight in time is determined from (A17) 
and (A13) where £t is the separation in hours between 
observation time and analysis time, the inherent temporal 
filtering of the analysis depends on the values of B t and 
Y £. It is desirable to compute B t from time correlations of 



different data types, however very little literature exists 
on this type of study. Barnes (1973) used the phase speed of 
the meteorological systems to compute values for that are 
consistent with the spatial parameters. Assuming a phase 
speed, C, of the order of 20 kt, then 8^. can be determined 
from B^/C, which is about 20 hours. For the update interval 
used, however, the temporal variations of the weights is 
only about 10%. 

The amount of horizontal and temporal filtering 
inherent in the analysis is governed by y^ and y t , 
respectively. The observational accuracy may be a factor in 
the determination of these values, as observations 
containing large random error tend to produce fictitious 
short waves which should be filtered more than the longer 
waves. The density of reports also becomes a factor, since 
aliasing may result when insufficient reports are used to 
describe a wave (Stephens, 1972). The filtering parameters 
selected (y H t = 0.3) give a response for the four-grid 
increment wave of only 25% with a fairly steep rise to 30% 
for the eight-grid increment waves. 

The error checking procedures used in analysis schemes 
may account for sizeable differences between various 
techniques. Even the most sophisticated schemes are not 
immune to large errors in the analysis that are caused by 
improper rejection of the observations. The most difficult 



problem is avoiding the rejection of observations that 
result from a poor forecast. This is most likely to happen 
if the forecast is used to check the data, rather than the 
more desirable "buddy check" method, in which observations 
are compared. In the design of the error rejection 
procedure, an attempt was made to retain as many 
observations as possible, even at the expense of accepting 
data with errors. These erroneous data cause small scale 
effects that tend to diminish during the balancing 
procedure . 

Three separate procedures are used to detect erroneous 
data. First, the radiosonde data are subjected to a 
vertical consistency check that requires the lapse rate be 
stable. Data is corrected through interpolation from 
adjacent levels with good data whenever possible. Secondly, 
gross errors are removed by rejecting the data that disagree 
with the forecast by more than five standard deviations of 
the expected error at that level. Finally, the remaining 
observations are used to perform a single pass, two- 
dimensional analysis. Each observation is then checked by 
first removing its effect from the analysis. This is done 
by determining the effect the observation in question has on 
the analysis at the closest grid point. 
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where is the weight given f r at this location. The 
analysis at this point excluding the observation being 
tested becomes 

Fq = F - Afp. (A28) 

The observation is rejected whenever disagreement between F c 
and fp exceeds three standard deviations of the forecast 
error and there exists the equivalent of two other observa- 
tions within one grid interval of the observation in 
question. This procedure rejects less than 0.2% of the 
total observations. An advantage of this method is that it 
utilizes computer code that is used to perform the analysis, 
which decreases the size and complexity of the computer 
program . 

The successive corrections procedure described above is 
both simple and very fast compared to other methods, such as 
the multivariate optimum interpolation method (Schlatter, 
1975). Furthermore, in experiments performed by Otto- 
Bliesner et al. (1977), the more sophisticated and time 
consuming methods did not produce significantly better 
results. Unfortunately, the successive corrections method 
is univariate, which means that no attempt is made to 
constrain the corrections to be consistent dynamically. 

This makes the balancing component of the assimilation 
system critical if the analyses are to be optimally combined 
to produce meteorologically consistent updates. 
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APPENDIX B 



LINEARIZED HYDROSTATIC, THERMODYNAMIC 
AND CONTINUITY EQUATIONS 



The vertical grid (see Fig. Bl ) has all variables except 
pressure are defined at odd levels and therefore 

4> k ~ ^k + 2 - ( "p^k + 2~ P k^ 0 k+l Tor k — 1 , 3 , . . . , k - 2 (Bl) 

and 

*k = *s + c p (p s‘ p k )9 k • < 82 ) 

Interpolation formulas 



1 1 



p 1 p 1 

p k+l ' H k -1 



< 1+k P k+] - P 



k - 1 



( B 3 ) 



1 1 
rr tt f 



p i p i "*’ k: 

p s H k-1 

P s - P k-1 



(B4) 



and 



e k+l = A k+1 0 k + B k + 1 0 k+2 



( B5 ) 



are used to produce energetically consistent equations where 
Pk 1 a k (Ps-Po*’ 



p k < 

p k - 



( B6 ) 
(B?) 



142 



A 



k+1 



B 



k+1 



P k+1 ~ P K. 
P k+2 " P k 

P k+2 " P k+1 
P k+2 " P k 



for k=l ,3 . . . k-4 



CBS) 



A 



K-i 



P K P K-2 

A 

P K-1 



K-2 



P K _ P k-2 



(B9) 



B 



K-l 




P K P K-2 




(BIO) 



T is temperature, p is pressure and C is the specific heat 

i 

at constant pressure. 



----------------- T l ^ j V L 

0 L+1 ’ P L+1 

T L + 1 ’ *L + 1 ’ V L+ 1 

FC ’ P SFC 

Fig. B 1 . Grid configuration. 
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The geopotential at each level is computed through 
integration, or 

*K = *S + C p (P S - P K> 9 K C BU > 

and 



♦k - + n f k C p (P n+2 - p n ) <W n ) + 



K, 



( B 1 2 ) 



2 C ( P - P ,) (B ,e ) 
n=k+2 p n n ' 2 n_1 n 



The primed sum indicates increments of 2. Now we can write 
the hydrostatic equation as 



K, 

*k " ^s n l } C kn T n 



where 



'kn 



0 




n<k 


W 2 


' p n^ A n+l /p n 


n = k 


c (p xo 
p n + 2 


- p n> A n + l/ p n 


+ 




C p (P n - P n 


-2> B n- 



or in matrix form 



? " ts = - I 



( B 1 3 ) 



n>k 



( B 1 4 ) 



The finite difference form of the thermodynamic 
equation (Eq. 299 in Arakawa and Lamb 1977) in orthogonal 
curvilinear coordinates is 
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[« e (F 



T s ) + 



V s + S7 p ij 5 > 9 ‘)?j 



r\ 

= ^-[(«c)|| + ip: v + (™>v + ” 9] ij ’ 

P (B15) 

where n = n , F = uu^p- , G = ttv^|- , S = na , the overbar 

is a linear average in the direction of the variable 
indicated, and <$ is a difference taken in the direction of 

A 

the subscript. 

To linearize, first subtract 




[ 5 1 n + 5 £ F + 



5 G + 
n 




0 , 



which gives 



ns 



t 





p Ws i]J a 





A a 



( B 1 6 ) 



( B 1 7 ) 




/ \ k 8 n _ 

(lTaa) ij n ' 



Qr) 



T'ij 



o r 



5 T k + rJlil p - T 1 

t i j L P,_ n K k ‘k J + 



k + 1 



k + 1 



ia. 



rliiii P 

lp k-l k 



- TJ 



k-1 



k J Aa i 



(7Taa) k 3ln W _ n , 

C p 3t " Vij 



( B 1 8 ) 



where is the rest-state temperature. 

Substituting the linearized form of the continuity 
equation (Eqs. 166-167 in Arakawa and Lamb 1972) give 

°k+l = - J; - V) n 40 n + J k + 1 jj 7 - ( V n )4 % * Q + 

( B 1 9 ) 
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k-1 



K 



k i = - s' (v-V) A a + a. E' ( V • V) Aa + Q 
k i ncl n n K 1 n =i n n 



K 

5 . 1 n it = - E'(v*V) Aa + Q 

t n =l n n p 



and 



T. dP, 

( ” oc,) k = C p (’’dT - ) 



gi ves 



+ n + 



k-2 K 

- 1) Z'(7-V) Aa + a. , E' (7*V) Aa 

_ n = l n n _ n= k + 2 n 

A a,. 



k-2 



- C3 



(o k-l ( 7 - V) n 4t, n + °k-l C 7 * v ) n ao n 



Aa, 



dP k K 



+ 3T-> n E,' (, ' v) n aa n ' « 



T ij 



In matrix form 



4 t T ,'j + = ( 9t>ij 



where 



r 



- n W+l" 1 ' n W-l" 1 * T k , dP k, . 

{C:l + da k - da k ' + (, 3rl* 4 »n " <k 



T kn = {C] + 



K+r 1 ) 

Aa ,, 



( °t+i ) 



a k-l 'k ^ P k 

- □_ + F 1 (^)>Aa 



Aa i 



p k 



n 



n = k 



r n k+1 > „ , a k-1 > , T k / dP k \ n . 

{[] + “aTT " ) + P7 17 d^r } } Aa n " >k 



( B20 ) 
( B21 ) 

(B22) 

n 



(B23) 
( B24 ) 



( B 2 5 ) 



The formula 
/» 



dP k _ f 

"JIT 



A k + 1 ^ a k + l (P k+2‘ P k^ + B k-1 ^°k-l ^ p k' P k-2^ 



( B26 ) 



k< K 



Aa i 



V p k 



for k=K 



is that derived by Arakawa from the interpolation formulas 
(B3) - ( B 1 0) . 

Note that the continuity equation may also be written in 
matrix form by defining 

rA0,“i 



1 



Aa- 



n = 



( B27 ) 



LAa„J 



so that 



5.1nir = -n6 + Q. 



( B28 ) 
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APPENDIX C 



EMPIRICAL ORTHOGONAL FUNCTIONS 



Holstrom (1963) showed that geopotent i al , <f>(x,y,p), 

may be accurately described in terms of a series in ortho- 
normal basis functions, 4> k (p), derived so that a partial 
sum. 



n 

(j> n ( x »y*P) = 2 a k (x,y)$ k (p) (Cl) 

k= 1 

where 



a k( x »y) = 1 <fr (x,y,p)i k (p)dp, (C2) 

h 

produces an optimum fit for all choices of n. In these 
equations, the horizontal coordinates are x and y and the 
vertical coordinate is p. The values p^ and P 2 are the 
vertical boundaries of the domain and the areal mean value 
of p at each level has been removed. 

Obukhov's (1960) method, which is computationally 
easier to use than Holstrom's (1963), uses the auto- 
covariance as a characteri Stic measure for determining the 
empirical orthogonal function, i.e., 

B ( P , p ' ) = Mx,y,p)*(x,y,p ' ) , (C3) 

where the overbar operator designates a horizontal average. 
This function describes the variance of p when p is equal to 
p', and it describes the covariance otherwise. The 
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redundancy of the ^-profiles is identifiable by the size of 
the covariance terms. Consequently, the accuracy of <f> n in 
(Cl) depends on the covariance magnitudes. For example, if $ 
is random, then (Cl) would not converge rapidly. Holstrom 
(1963), however, has observed considerable redundancy in the 
atmosphere for geopotential. 

In Obukhov's (1960) method, the eigenfunctions of the 
integral operator, 

P2 B(p,p ' )<f> (p ' ) dp ' = u k 4> (p) , (C4) 

1 

produce the optimum choice for basis functions, <j> k (p), which 
produce least mean square error for all values of n in (Cl). 
The eigenvalues of ( C 4 ) , y k , measure the variance that their 
associated eigenfunctions represent. Therefore, the order- 
ing of the functions is made so that p k is in descending 
order. Also, the averaged normalized variance of <1 > that is 
represented by 4> n is simply 




a 



n 




( C 5 ) 



To extend Obukhov's (1960) procedure to finite differ- 
ences, the independent function is represented by 

= i>( lAx, jAy.kAp) , 



i » j » k 



(C6) 



where ax, Ay, and Ap are the grid spacings and i,j,k the 
grid location. Then the autocovariance becomes 



B 



z ,m 



M N 

s z 



*1,0,* *i ,j ,m 



i=l 0=1 

The integral operator (C4) now has the form 
K 

Z B„ d> = y , d> . , 

£,nrm ^k y £ ’ 

m= z 



( C 7 ) 



where z and m are vertical indexes. The optimum basis 
functions are the eigenvectors of B z m arranged so that 
corresponding eigenvalues are in descending order. There 
are K modes or eigenvectors in this system, which correspond 
to the size of the square array Bji m . However, the 
motivation of this approach is to allow a partial sum to be 
used that maintains most of the accuracy of the original 
function .• p . 
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